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I.  Introduction 


This  research  focuses  on  developing  enhanced  contrast  thermal  acoustic  imaging  (TAI) 
technology  for  the  detection  of  breast  cancer  by  combining  amplitude-modulated  (AM) 
electromagnetic  (EM)  field  excitation,  re  sonant  acoustic  scattering,  and  advanced 
signal  processing  techniques.  EM- induced  TAI  com  bines  the  m  erits  of  both  EM 
stimulation  and  ultrasound  im  aging,  while  overcom  ing  their  respective  lim  itations. 
EM  imaging  provides  excellent  contrast  between  cancerous  and  normal  breast  tissue, 
but  the  long  wavelengths  provide  poor  sp  atial  resolution.  Conventional  ultrasound 
imaging  possesses  very  fine  m  illimeter-range  sp  atial  reso  lution  but  poor  soft  tissue 
contrast.  While  EM-induced  TAI  possesses  great  promise,  the  thermal  acoustic  signals 
tend  to  be  weak.  However  ,  whe  n  the  turn  or  is  exc  ited  into  reson  ance  via  E  M 
stimulation,  the  effective  acoustic  scattering  cross-section  may  increase  by  a  factor  in 
excess  of  100  based  on  predic  tions  for  microsphere -based  ultrasound  contrast  agents. 
Such  an  in  crease  wou  Id  truly  be  re  volutionary,  m  aking  the  EM-induced  T  AI 
technology  a  very  promising  candidate  for  routine  breast  cancer  screening.  To  induce 
the  resonant  response  from  the  tumor,  we  consider  various  approaches  including,  for 
example,  AM  continuous  wave  (CW  )  EM  stim  ulation,  where  the  m  odulation 
frequency  range  contains  the  predicted  re  sonant  frequencies  for  a  distribution  of 
tumor  sizes  and  contras  t  ratios.  The  carrier  frequency  of  the  EM  stimulation  can  be 
fixed  and  chosen  for  the  best  penetrati  on  and  heat  absorption.  The  im  age  formation 
methods  in  the  existing  T  AI  system  s  are  predom  inantly  data-independent 
delay-and-sum  (with  or  without  w  eighting)  type  of  approaches.  These  approaches 
tend  to  have  poor  resolution  (relative  to  the  best  possible  resolution  a  transducer  array 
can  offer)  and  high  sidelobe  problem  s,  especially  when  the  transducer  array  is  not 
composed  of  unifor  mly  and  linearly  spaced  transducers,  which  is  th  e  case  for  the 
existing  TAI  systems.  We  devise  adaptive  image  formation  algorithms  to  achieve  high 
resolution  and  excellent  interference  and  noise  suppression  capability. 
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II.  Body 


II.  1  Theory 


Previous  work  in  the  area  of  thermo-acous  tic  imaging  all  utilized  hig  h  peak  power , 
short  pulse  excitation  [1-4],  In  essence  th  ese  approaches  are  time  domain  based,  and 
capture  the  electro-acoustic  impulse  response  of  the  phantom  system  to  short,  intense 
EM  ilium  ination  with  a  prescr  ibed  (high)  energy  density .  T  hey  genera  lly  require  a 
prohibitively  expensive  power  am  plifier  along  with  broadband  m  icrowave 
components  and  a  high-speed  data  acquisiti  on  system.  The  present  study  approaches 
the  problem  from  the  frequency  dom  ain,  an  d  seeks  the  sam  e  inf  ormation  as  the 
time-domain  approaches  but  with  using  lower  power,  narrow  band  excitation  to  obtain 
the  steady-state  response.  It  is  anticipated  that  due  to  the  thermo-acoustic  resonance  of 
phantom,  similar  imaging  information  can  be  obtained  but  with  a  significant  reduction 
in  excitation  power. 

The  theoretical  analysis  that  forms  the  foundation  for  Thermo-Acoustic  Imaging  (TAI) 
is  based  on  Diebold’  s  theory  publishe  d  in  1988  [5].  There  he  provides  the 
theoretical  analys  is  for  pressure  wave  generation  by  exciting  droplets  with  a 
modulated  laser  pulse.  Using  Diebold’s  approach,  it  can  be  shown  that  the  steady-state 
pressure  response  of  the  phantom  to  an  AM  EM  wave  will  be  of  the  form: 


Pf 


jj3crE2csa 


(sing -geos  g)-^- 
9 


vv 


A_ 

Pf 


S— ^ -  cos  q  +  y'^^sing 
g  pfcf 


~J<F 


(1) 


Where  the  phantom  has  the  properties: 

P  =  thermal  expansion  coefficient  [  1/K] 

a  =  electrical  conductivity  [S/  m] 

Cp  =  specific  heat  [J/kg*  K] 

E  =  electric  f  ield  intens  ity  inside  th  e  phantom  (as  sumed  to  be  unifor  m 
over  the  volume  of  the  phantom)  [Y/m] 

cs  =  speed  of  sound  in  phantom  [m/s] 

cf  =  speed  of  sound  in  surrounding  material  [m/s] 

ps  =  density  of  phantom  [kg  /  m3] 
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pf  =  density  of  surrounding  materials  [  kg  /  m3  ] 

a  =  radius  of  phantom  [m] 

Additionally  the  term  s  q,  f  ,  and  r  are  the  norm  alized  m  odulation  frequency  , 
normalized  retarded  time,  and  normalized  position  of  transducer,  respectively,  and  are 
defined  as: 


-V  c 
t  =—t 
a 


t  =  — 
a 


"/ 


coa  7 

q  =  —  =  ksa 
c 


r  =  ■ 


dimensionless  time 


retarded  time 


wave  vector 


distance  from  sphere  center 


From  this  analysis,  we  can  conclude  that  the  frequency  dom  ain  response 

characteristics  of  the  phantom  pressure  signal  are  determined  primarily  by  the  density 
ratio  and  sound  speed  ratio  betw  een  phantom  and  surrounding  materials.  The 
amplitude  of  pressure  response  prim  arily  relies  on  the  m  aterial  property  of  phantom 
(thermal  expansion  coef  ficient,  conductivity ,  sp  ecific  heat),  electric  field  intens  ity 
imposed  on  phantom,  and  relative  position  of  transducer. 

This  analys  is  constitu  tes  the  foundation  of  our  num  erical  sim  ulation.  Until  very 
recently  we  have  not  had  a  high  degree  of  confidence  in  two  critical  param  eters, 
namely  the  density  of  p  hantom  and  the  elect  ric  field  intensity  distribution  which  can 
be  obtained  from  EM  field  m  easurements  in  the  tank.  For  other  m  aterial  constants, 

specifically  ft ,  cr ,  C  ,  and  cs ,  w  e  have  re  lied  on  in  telligent  estim  ates  and  are 

currently  in  the  pro  cess  of  designing  experiments  to  accurately  measure  and  confirm 
the  values  used.  Figure  1  shows  the  expected  phantom  pressure  response  as  a  function 
of  modulation  frequency  based  on  the  assumed  material  constants. 
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Figure  1 :  Phantom  pressure  response  as  a  function  of  modulation  frequency. 


Note  that  the  lar  gest  r  esonant  pea  ks  o  ccur  at  frequencies  which  are  outside  the 
bandwidth  of  the  transducer  used  (center  frequency  1MHz,  frequency  band  0.6  MHz, 
Panametrics  NDT  ,  Model:  Y303).  Hence,  in  addition  to  better  determ  ining  the 
material  properties  of  the  phantom  ,  a  smaller  phantom  with  hi  gher  resonant 
frequencies  is  currently  being  designed  and  fabricated. 


Although  Diebold’ s  analysis  was  carried  out  in  frequency  dom  ain,  the  experim  ents 
were  completed  using  very  short  laser  pul  ses  that  are  well-approxim  ated  as  a  Dirac 
delta  function.  Consequently  the  sinusoidal  steady-state  analysis  approach  followed  in 
the  present  research  ef  fort  para  lleled  th  e  D  iebold  solu  tion.  There  is,  how  ever,  a 
significant  difference  in  the  two  approaches .  It  was  pointed  out  by  Lihong  Y.  Wang  in 
2000  [6]  that  if  the  thermal  confinement  condition  does  not  apply,  the  heat  conduction 
effects  should  be  taken  into  consideration  in  calculating  the  pre  ssure  wave  generated 
by  electrom  agnetic  ilium  ination.  In  our  case,  a  Continuous  W  ave  (CW)  m  odulated 
microwave  source  was  used  to  heat  the  pha  ntom  and  was  applied  over  a  m  uch  longer 
time  scale  as  com  pared  with  im  pulsive  excitation.  Conseque  ntly  heat  conduction 
effects  must  necessarily  be  included  in  the  analysis. 


A  sim  pie  num  erical  simulation  demonstrat  es  the  ef  fects  due  to  heat  conduction 
between  phantom  and  the  surrounding  water.  Consider  the  simple  model  illustrated  in 
Figure  2. 
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P  =  Phantom 
W  =  Water 

Smaller  sphere  radius  a 
=  2.4  mm 

Larger  sphere  radius  b 
=  60  mm 


Figure  2:  Model  of  the  phantom  and  the  surrounding  water. 

The  phantom  volume  is  heated  by  heating  function  G(t,  r),  which  is  uniform  over  the 
volume  of  the  phantom  and  results  from  the  elec  tromagnetic  illu  mination.  The 
boundary  condition  between  phantom  and  water  is  that  th  e  temperature  and  heat  flux 
in  phantom  and  the  water  be  continuous  at  the  phantom  -water  interface.  The  water 
volume  is  assumed  to  be  very  large  hence  the  outside  boundary  of  water  (the  water-air 
interface)  is  assumed  to  be  perfectly  insulated.  The  governing  equations  are: 


1  d  (J2  dTp(r,t)  ap  dTp(r,t ) 


ap— : - (r 

P  r2  dr 


dr 


~)  +  ~r~Sp{r  ,t)  =  - 


1  d  2dTw(r,i)  dTw(r,i) 
w  r 2  dr  dr  dt 

Tp(0,t)  =  finite 

Tp(a,t)  =  Tw(a,t) 


'IV 


dTp(r,t)  | 

dr 

dTw(r,t) 
dr 


-  k 

r=a  ^ W 


dr 


r=b 


=  0 


dt 


0  <  r  <  a 


a  <r  <b 


(2) 


In  Equation  (2),  ap ,  kp  are  th  ermal  characteristics  of  phantom  ,  aw  ,kw  are  th  ermal 
characteristics  of  wat  er,  Tp  is  temperature  distribution  in  phantom  ,  Tw  is 

temperature  distribution  in  water,  and  gP(r,t )  is  the  heating  function  resulting  fro  m 

the  electromagnetic  illumination.  In  order  to  solve  thes  e  equations,  the  appropriate 
Green’s  function  has  been  found  and  utilized.  The  solution  is  in  the  form  of  series 
expansion  and  is  given  by  (note  that  500  terms  are  used  in  the  expansion): 
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Since  the  tem  perature  distribution  in  the  phantom  and  water  is  known,  the 
temperature  distribution  can  be  determ  ined  and  is  shown  in  Figure  3.  In  Figure  3  the 
solid  line  represents  the  radial  temperature  distribution  from  the  center  of  the  phantom 
to  its  boundary  (5  mm ),  while  the  dashed  li  ne  represents  the  te  mperature  distribution 
when  only  the  phantom  being  heated  and  not  the  surrounding  water  .  The  dif  ferent 
colors  illustrate  how  the  temperature  distribution  evolves  as  a  function  of  time.  From 
Figure  3  it  is  seen  th  at  the  tem  perature  in  the  phantom  is  being  reduced  at  the 
boundary  due  to  heat  conduction  to  surrounding  water.  This  effect  acts  as  a  negative 
factor  for  steady-state  electromagnetic  illumination. 


Comparison  of  temperature  distribution  in  two  situations 
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Figure  3:  Temperature  distribution. 

Generally  speaking,  the  temperature  in  phantom  will  present  a  linear  increasing  trend, 
together  with  a  sm  all  sinusoidal  oscillation  upon  it,  caused  by  th  e  modulating  signal. 
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The  governing  equation  describing  the  acoustic  pressure  wave  is  the  wave  equation, 
which  can  be  expressed  as: 


1  d2ps  _  -p  dH 
]  ^  “  Cp  dt 


In  (5),  ps  is  the  sound  pressure,  cs  is  the  sound  speed,  P  is  the  th  ermal  expansion 

coefficient,  Cp  is  the  h  eat  capacity,  and  H  is  the  h  eating  function  of  the  phantom.  In 

solving  this  equation,  the  heat  conduction  issue  discu  ssed  above  should  be  taken  into 
consideration,  namely  that  the  heating  function  should  be  expressed  as: 

H(r,  t)  =  ^4(1  +  cos  cot)  -  H  \r,  t)  (6) 


The  first  item  on  the  right  side  of  Equation  (6)  represents  the  electromagnetic  heating, 
while  the  second  item  represents  the  heat  conduction  into  water.  Including  the  thermal 
conduction  considerations  discussed  Equation  (7)  must  satisfy: 


and  hence 


H(r,t)  =  pCp 


dTjrJ) 

dt 


(7a) 


dH{r,t)  _  a2T(r,Q 
dt  P  p  dt 2 


(7b) 


At  this  point,  it  is  assumed  that  the  temperature  distribution  is  a  linear  function  of  time 
with  a  small  superimposed  sinusoidal  oscillation.  This  may  be  expressed  as: 

T  (r,  t )  =  A(r)t  +  AT(a>)  cos  (cot)  (8a) 


giving 


dH(f,  t)  _  pc  \T{co)(o2  cos  {cot)  (8b) 

dt 

Applying  the  Diebold  solution  one  obtains: 

dH(rJ)  =  aI()o)^-JO,t  =>  al(]  =  pCAT{o))o)  (9) 
dt 

Equation  (9)  can  now  be  substituted  into  the  Diebold  solution  resulting  in  the  acoustic 
pressure  wave  expressed  as  the  result  of  the  sinusoidally  increas  ing  tern  perature 
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component,  namely, 


Pf 


joj{3cspsaAT(co) 


(sing-gcosg)-^- 

q 


vv 


1_A 

Pf, 


Sm  —  -  cos  q  +  j  ^',C'  sin  q 
q  pfcf 


(10) 


The  temperature  increase  AT  can  be  estim  ated  from  the  heat  conduction  analysis 
above. 

It  is  difficult  to  analytically  decouple  the  pressure  wave  generation  excited  by  the  CW 
excitation  and  heat  conduction  ef  fects.  It  is  apparent  from  the  above  analysis  that  for 
the  case  of  CW  excitation,  the  he  at  c  onduction  ef  fects  will  grea  tly  attenua  te  the 
electromagnetic  ilium  ination.  The  only  w  ay  to  m  itigate  this  a  ttenuation  is  to  u  se 
significantly  lar  ger  CW  excita  tion.  Short  of  doing  this,  the  signa  1  is  g  enerally  too 
weak  to  be  detected. 


II.2  Phantom  Development 

(a)  Phantom  Modeling 

A  tissue  phantom  for  therm  al  acoustic  im  aging  should  m  atch  the  living  tissue’  s 
electrical  as  well  as  aco  ustic  properties.  The  preliminary  phantom  development  is 
focused  on  m  atching  dielectric  properties,  specifically  permittivity  and  conductivity. 
According  to  Duck  [7],  the  relative  permittivity  of  malignant  breast  tissue  at  a  range 
encompassing  434  MHz  ranges  between  36-56,  and  the  conductiv  ity  ranges  between 
0.35-0.8  S/m.  The  phantom  was  m  ade  from  TX-151  powder  (f  rom  Oil  Research 
Center),  tap  water  ,  cane  sugar  and  pot  assium  chloride,  accordin  g  to  m  ethods 
developed  at  The  McKnight  Brain  Institute  at  the  University  of  Florida  by  B  eck  [8]. 
The  dielectric  properties  are  controlled  by  varying  the  concentration  of  the  cane  sugar 
and  potassium  chloride.  An  HP85070B  coaxial  probe  and  HP8752C  network 
analyzer  measured  the  dielectric  properties.  A  preliminary  sample  was  created  using 
the  experien  ce  of  the  staf  f  of  the  McKnight  Brain  Institute  in  creating  brain  tissue 
models  that  have  diele  ctric  prope  rties  s  imilar  to  that  of  m  alignant  tissue.  A  n 
iterative  process  of  m  easuring  the  diele  ctric  constan  ts  and  m  anipulating  the  cane 
sugar  and  potassium  chloride  concentration  was  conducted  until  a  suita  ble  phantom 
was  created.  The  tissue  phantom  created  is  used  in  the  excitation  system.  A  diagram 
of  the  experim  ent  is  shown  in  Figure  4.  The  cylinder  height  is  10.3  inches  and  the 
inner  diameter  is  7.25  inches.  The  transducer  used  in  the  experim  ent  is  a  0.5  inch 
diameter  piezoelectric  immersion  transducer  from  Olympus  NDT  (model  V303).  It 
has  a  center  frequency  of  0. 91  MHz  and  60.275%,  -6  dB  bandwidth.  At  1  MHz,  the 
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half  angle  beam  width  is  approxim  ately3.5°.  R  eciprocity  calibra  tions  w  ill  be 
performed  in  de-ionized  water  .  The  wa  veguide  used  is  from  Penn  Engineering 
Components  (model  WR  187)  with  a  frequency  range  of  3.95-5.85  GHz. 

0 


E-field  probs 


Figure  4:  Phantom  experimental  setup. 


(b)  Initial  Phantom  Resonance  Experiments 

Experimentation  designed  to  verify  acous  tic  resonance  of  the  turn  or  phantom  by 
electromagnetic  excitation,  with  the  experiment  setup  shown  schematically  in  Figure 
1 ,  consists  of  an  electrom  agnetic  wavegui  de,  two  ultrasonic  transducers  (one  for 
transmission  ant  the  other  for  sensing),  a  nd  a  Plexiglas  cy  linder  that  contain  s  the 
de-ionized  water  and  tumor  phantom.  The  cylinder  height  is  10.3  inches  and  the  inner 
diameter  is  7.25  inches.  The  transducer  used  in  the  experiment  is  a  0.5  inch  diam  eter 
piezoelectric  immersion  transducer  from  Olympus  NDT  (model  Y303).  It  has  a  center 
frequency  of  0.91  MHz  and  60.275%,  -6  dB  bandwidth.  At  1  MHz,  the  half  angle 
beam  width  is  approxim  ately  3.5°.  Reci  procity  calibrations  are  perform  edin 
de-ionized  water.  The  waveguide  used  is  from  Penn  Engineering  Components  (model 
WR  187)  with  a  frequency  range  of  3.95-  5.85  GHz.  In  the  experim  ents,  a  1  cm 
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diameter  sample  of  the  phantom  is  suspended  in  the  center  of  the  cylinder,  submerged 
in  de-ion  ized  w  ater,  and  aligned  w  ith  the  ultra  sonic  tr  ansducers.  The  transm  itting 
transducer  is  mounted  to  the  side  of  the  cylinder  with  the  trans  ducer  face  exposed  to 
the  de-ionized  water.  The  rece  iving  transducer  is  similarly  aligned  and  its  axis  m  akes 
a  180-degree  with  respect  to  the  first  transducer.  The  electrom  agnetic  waveguide 
provides  excitation  from  the  bottom  of  the  cylinder.  Electric  field  probes  are  inserted 
through  a  port  in  the  top  of  the  cylinder  and  are  used  to  measure  total  field  intensity 
along  the  cylinder  axis.  Two  probes  are  need  ed:  one  to  m  easure  the  radial  field 
intensity  and  another  to  determine  the  axia  1  field  strength.  The  ra  dial  field  probe  is 
shown  in  Figure  8.  The  goal  of  the  prelim  inary  setup  is  to  achieve  a  uniform  field 
within  a  3  cm  radius  around  the  phantom  .  The  frequency  of  the  434  MHz  radio 
frequency  signal  can  be  amplitude  modulated  over  a  range  on  the  order  of  400  kHz  to 
4  MHz,  which  is  th  e  expected  ran  ge  the  resonant  frequency  of  the  phantom  .  The 
acoustic  signal  from  the  phantom  is  m  easured  by  the  ultrasonic  transducer  and  is 
recorded  by  the  NI  data  acquisition  system. 


Figure5:  Test  setup  for  determining  the  resonant  frequency  of  the  tumor  phantom. 

Initially,  the  turn  or’s  resonant  frequenc  y  is  established  acoustically.  Referring  to 
Figure  5,  a  swept  sinusoidal  signal  generated  emanating  from  Port  1  of  an  RF  network 
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analyzer  acts  as  a  source  which  drives  one  acoustic  tran  sducer  while  the  secon  d 
acoustic  transducer,  used  to  receive  the  acoustic  signal,  drives  port  2.  W ith  the  tumor 
absent  from  the  wate  r- filled  tank,  the  receiv  ed  signa  1  c  onsists  of  the  acoustic 
frequency  characteristics  of  the  tank  itself,  and  this  frequency  response  is  used  as  a 
reference  to  calibrate  the  network  analyzer.  The  same  sweep  is  then  perform  ed  with 
the  tumor  phantom  present  and  the  response  no  w  contains  the  resonant  characteristics 
of  the  phantom .  A  typical  frequency  respons  e  is  shown  in  Figure  6  where  it  is  seen 
that  the  resonant  frequency  occurs  between  around  1  MHz. 

For  TAI  sys  terns  the  acoustic  pres  sure  wave  is  generated  as  the  tumor  is  heated  via 
application  of  EM  energy.  Having  establishe  d  the  resonant  frequency  of  the  phantom 
using  the  a  coustic  tech  nique  descr  ibed,  a  sig  nal  w  ith  th  is  f  requency  am  plitude 
modulated  at  434  MHz  EM  carrier  signal  is  then  used  to  drive  the  waveguide  located 
at  the  botto  m  of  the  tank  via  an  RF  power  amplifier.  The  acoustic  transducer  is  the  n 
used  to  sense  the  acoustic  pressure  wave  generated  by  the  TAI  process. 
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Figure  6:  A  typical  frequency  response  of  acoustic  resonance. 


(c)  Phantom  Aging  Characteristics 

While  carry  ing  out  these  m  easurements,  it  becam  e  apparent  the  large  num  her  of 
variables  in  volved  in  c  haracterizing  the  resonance  of  the  phantom  ,  one  being  the 
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effect  of  aging  of  the  phantom  ove  r  tim  e.  Subsequently,  som  e  ti  me  was  focused  to 
characterize  the  phantom  aging.  One  of  the  aspects  explored  was  the  v  ariation  of  the 
resonant  frequency  over  time.  Figure  7  shows  results  based  on  the  experimental  setup 
described  in  the  previous  section  and  s  hown  in  Figure  5.  W  hile  the  sam  pie  size 
consists  of  only  3  phantoms,  the  general  trend  seems  to  be  that  the  resonant  frequency 
initially  increases  before  it  decreases  over  time.  More  specimens  need  to  be  studied  to 
confirm  this  trend.  In  Figure  7,  phantom  c  exhibited  no  apparent  resonance  in  week 
one.  This  was  attributed  to  phantom  pi  acement  in  the  cylinder.  The  setup  was 
consequently  altered  in  a  way  that  would  ins  ure  repeatable  phanto  m  placem  ent. 
Phantom  aging  is  important  to  quantify  in  order  to  minimize  experimental  error  over 
time.  Other  properties  are  currently  being  examined  including  density  and  speed  of 
sound. 


Resonant  frequency  as  a  function  of  time 
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Figure  7:  Resonant  frequency  of  phantom  over  time.  Note  no  resonance  was  found  for 

phantom  c  in  the  first  week. 


(d)  Electric  Field  Measurements 

Accurate  modeling  of  both  the  acoustical  and  electrical  behavior  of  the  system  forms 
an  essential  com  ponent  of  this  research  e  ffort.  Of  special  concer  n  is  the  ability  to 
model  ele  ctrical-to-acoustical  tran  sduction  asp  ect  of  this  system  .  The  acous  tic 
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pressure  wave  to  be  measured  is  generated  when  the  modulated  electrical  signal  heats 
the  phantom  causing  it  to  mechanically  undulate.  This  heating  occurs  as  a  result  of  the 
electric  field  intensity  in  the  location  of  the  conducting  phantom  .  As  a  consequence 
not  only  must  the  electrical  fields  p  resent  in  this  system  be  accu  rately  modeled,  bu  t 
these  fields  must  also  be  verified  so  that  the  proper  heating  parameters  may  be  used  as 
input  to  the  acoustic  m  odel.  Initially  an  off-the-shelf,  low-pr  ofile,  high-resolution 
probe  that  measures  the  electric  potential  (Carsten  Associates  Model  E-601)  was  used. 
Once  the  potential  is  kn  own  a  differencing  procedure  can  be  used  to  determ  ine  the 
field.  It  was  found,  however,  that  this  method  proved  inadequate  for  the  present 
system.  In  particular,  the  common  m  ode  signal  m  easured  far  outweighed  the 
difference  signal  and  an  accurate  estimate  of  the  electric  field  could  not  be  determined. 
Further,  the  metal  shaft  of  the  probe,  though  small  in  diameter,  was  distorting  the  field 
in  the  tank.  To  remedy  these,  two  custom  probes,  one  used  to  measure  the  axial  field, 
another  used  to  measure  the  radial  field  were  designed.  Since  the  new  probes  measure 
the  electric  field  directly  no  differencing  is  required.  Furthermore  the  new  probes  are 
design  with  a  shaft  that  is  impedance-matched  to  the  im  pedance  of  de-ionized  water 
thus  making  the  probe  transparent  except  of  course  (unavoidably)  for  the  imm  ediate 
region  being  measured.  To  date  the  radial  field  probe  has  been  obtained  and  is  shown 
in  Figure  8.  Shown  also  in  Figure  8  is  a  Tr  ansverse  Electro-Magnetic  (TEM)  test  cell 
used  for  probe  calibration.  One  end  of  the  TEM  cell  is  excited  with  a  434  MHz  signal 
while  the  other  end  is  terminated  in  a  50-ohm  load.  Since  the  dimensions  of  the  cell 
are  precisely  known,  the  electri  c  field  in  the  cell’s  center  is  also  known.  This  center 
portion  of  the  TEM  cell  is  filled  with  de-ionized  water  in  which  the  probe  is  inserted, 
and  the  res  ulting  vo  ltage  is  m  easured  using  a  50-ohm  oscillos  cope.  O  nee  the  axia  1 
filed  probe  is  com  pleted,  an  accurate  field  map  of  the  tank  along  the  ax  is  where  the 
phantom  will  be  p  laced  will  be  obtained.  Once  again  this  data  will  be  used  to  verify 
the  electromagnetic  model  and  link  the  electrical  and  acoustic  models. 
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(e)  Electrical-to-Acoustical  Transduction 


To  date  two  problems  have  been  encountered  with  regard  to  accurately  measuring  the 
RTAI  effec  t  and  are  currently  being  rem  edied.  The  first  regards  drift  of  t  he 

acousticacoustic  resonance  m  easurements.  Since  the  acou  stic  tran  sducers  poss  ess 
narrow  beamwidths,  accurate  placem  ent  of  the  phantom  is  required  f  or  a  stable 
resonance.  Secondly  the  acousti  c  resonance  varies  with  the  temperature  of  the  water. 
Initially  the  acoustic  and  electrical  porti  ons  of  the  experim  ental  procedure  where 
performed  consecutively.  It  has  been  de  termined  that  m  onitoring  the  acoustic 
resonance  while  the  electrical  m  easurements  are  taken  is  essential.  F  urthermore,  to 
date  the  electrical  (RF)  power  level  has  been  kept  low  with  a  total  drive  power  of  less 
than  four  watts.  This  requires  that  the  measurement  bandwidth  be  kept  small  so  as  to 
not  allow  excessive  noise  in  the  measurements.  A  new  RF  power  amplifier  capable  of 
producing  a  drive  power  of  up  to  100-watts  has  been  obtained.  This,  along  with 
accurate  p  lacement  of  the  phantom  ’s  location  will  a  llow  for  a  robust,  s  table 
measurement  of  the  RTAI  effect. 


II.3  Electromagnetic  Stimulation 

(a)  Tumor  Phantom 

For  TAI  systems  the  acoustic  p  ressure  wave  is  generated  as  the  turn  or  is  heated  via 
application  of  EM  ener  gy.  For  a  proof-of-  concept  leve  1  dem  onstration  of  the  T  AI 
process,  the  EM  excitation  is  acco  mplished  with  the  experim  ental  setup  shown  in 
Figure  9.  A  Plexiglas  tank  in  the  shape  of  a  right  circular  cylinder  (length  =  10.29 
in.,  radius  =  4  in)  is  filled  with  de-ionized  water  with  relative  permittivity  s/e0  =  81, 
and  conductivity  a  «  0.  The  quality  of  the  reconstruc  ted  image  is  direc  tly  related  to 
the  uniformity  of  the  EM  excitation  in  the  cavity.  For  the  cavity  dimensions  used  here, 
the  EM  f  ield  will  be  th  e  sum  of  several  cavity  modes,  and  a  unif  orm  EM  f  ield  can 
only  be  approximated  with  multiple  exciters.  For  the  present  case,  a  single  source  is 
used  to  exc  ite  the  cavity  so  a  uniform  field  is  not  exp  ected.  The  RF  excitation  is 
achieved  using  a  coax- to- waveguide  ad  apter  (WR-187)  which  has  an  operating 

frequency  range  of  3.95  -  9.85  GHz.  Though  th  e  frequency  of  the  exciting  signal  is 
434  MHz,  a  value  typical  for  m  edical  im  aging  system  s,  to  m  inimize  m  ismatch 
between  the  waveguide  exit  aperture  and  the  water-filled  tank,  the  waveguide  itself  is 
water-filled,  thus  reducing  the  wavelength  in  the  w  aveguide  to  3.0  in.,  andth  us 
operating  the  exciting  waveguide  within  its  recommended  range. 

The  tumor  is  sim  ulated  by  m  eans  of  a  m  embrane  filled  with  a  gel  that  p  rovides  the 
desired  permittivity  and  conductivity.  This  turn  or  will  be  s  uspended  from  the  top  of 
the  cavity  along  the  axis  of  th  e  cylinder  at  a  location  of  a  local  maximum  of  electric 
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field.  The  location  of  the  maximum  electric  field  is  determined  by  using  a  low-profile, 
high-resolution  E-field  probe  (C  arstens  Associates  Model  E-601).  The  probe 
impedance  is  50  Q  and  the  relative  s  ignal  intensity  is  measured  with  an  RF  spectrum 
analyzer.  The  probe  tip  is  tr  anslated  along  the  ax  is  of  the  cylinder  and  the  location(s) 
of  the  maximum  field  is  determined.  The  tumor  is  then  placed  at  this  location  with  the 
E-field  probe  removed. 


Figure  9:  Schematic  representation  of  the  experimental  setup  used  to  excite  and 
measure  the  cavity’s  axial  electric  field.  Also  shown  is  the  electric  field  obtained 
under  assumed  operating  conditions  using  a  FDTD  simulation. 

The  constru  ction  of  th  e  tes  t  se  t-up  as  desc  ribed  is  com  plete.  Th  e  R  F  signal  is 
generated  with  a  signal  s  ynthesizer  am  plified  with  a  power  am  plifier  capable  of 
providing  an  output  power  on  the  order  of  10  watts.  The  acoustic  signal  generated  by 
the  tumor  expansion  will  be  detected  with  a  pressure  sensor(s)  as  described  elsewhere 
in  this  report. 

Anticipated  ele  ctromagnetic  perf  ormance  has  been  m  odeled  by  m  eans  of  a  Finite 
Difference  T  ime  Do  main  ( FDTD)  si  mulation  program .  Si  nee  the  size  of  the  cavity 
used  along  with  the  fac  t  that  it  is  filled  with  a  h  igh  permittivity  material  suggests  that 
several  cavity  m  odes  will  exist,  with  the  total  f  ield  equaling  the  sum  of  these  m  odal 
fields.  Electromagnetic  modeling  along  with  due  consideration  of  the  cavity’s  acoustic 
properties  provided  the  guidelines  used  to  d  etermine  th  e  cavity  d  imensions.  For 
example,  the  specific  cylinder  leng  th  and  di  ameter  was  selected  in  part  by  ensu  ring 
that,  via  num  erical  simulations,  the  electric  field  was  not  at  a  spatial  null  along  the 
cylinder  axis  and  that  at  least  one  electr  ic  field  peak  occurred  along  this  axis.  For 
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simulation  purposes,  the  high  dielectric  permittivity  -  air  interface  was  modeled  by  a 
perfect  magnetic  conductor  (PMC)  at  all  cavity  boundaries.  The  relative  amplitude  of 
the  electric  field  is  shown  along  side  the  cavity’s  schematic  representation  in  Figure  9. 


(b)  Hemispherical  breast  model 

Use  of  the  RTAI  concept  as  a  real  im  aging  application  requires  extending  the 
approach  to  a  three-dimensional  system.  As  the  experimental  portion  of  present  phase 
of  the  research  effort  co  ntinues  so  does  the  modeling  for  the  next  phase  of  the  effort 
which  involves  three-dimensional  modeling,  simulation,  and  measurements. 


Excitation  voltage:  +v,(0- 


Figure  10:  Motivation  for  achieving  a  uniform  field  within  a  hemispherical  region 
using  small,  capacitively  coupled  patches. 

Figure  10  illustra  tes  co  nceptually  how  th  e  required  unif  orm  electric  field  for  a 
realistic  breast  phantom  can  be  achieved  in  a  simple  and  practical  m  anner.  Consider 
the  situation  shown  in  Figure  10(  a)  which  illustrates  a  diele  ctric  sphere  of  arbitrary 
permittivity  in  free  space.  It  is  well  known  in  electrom  agnetic  theory  tha  t  a 

sinusoidally  distributed  surface  charge  dis  tribution,  ps  =  p{)  cos  6 ,  w  ith  p()  as  a 

constant,  placed  on  the  surface  of  the  sphere  results  in  a  perfectly  un  iform  electric 
field  inside  the  dielectric  sphere  as  shown.  By  invoking  the  theory  of  images  a  Perfect 
Magnetic  Conductor  (PMC),  which  represents  a  good  model  for  the  chest  wall,  can  be 
placed  as  shown  without  altering  the  electric  field  in  the  sphere.  In  this  way  a  uniform 
electric  field  in  a  hemispherical  region  that  models  the  human  breast  can  be  obtained. 
Obtaining  such  a  surface  charge  distributio  n  is,  in  general  not  possible,  though  it  can 
be  approximated  in  a  pointwise  fashion  as  illustrated  in  Figure  10(b).  There  it  is  seen 
that  small  (subresonant)  patche  s  can  be  placed  o  n  the  breast  so  as  to  approximate  a 
uniform  electric  field  inside  the  breast.  These  (metal)  patches  would  be  coated  with  a 
non-conducting  material  and  hence  modeled  as  capacitively  coupled  as  shown.  The 
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concept  is  akin  to  the  approach  used  to  obtain  a  uniform  ,  circularly  polarized 
transverse  m  agnetic  field  from  discrete  wire  cage  stru  cture  in  Magn  etic  Resonance 
Imaging  (MRI)  applications.  Though  the  perf  ectly  uniform  field  would  be  obtained 
only  for  the  case  of  an  electros  tatic  charg  e  distribution,  it  is  not  unreasonable  to 
expect  a  quasi-static  solution  which  reasonably  approximates  the  ideal,  static  situation 
under  appropriate  conditions,  and  this  is  cu  rrently  being  explored.  To  illustrate  the 
basic  idea,  Figure  1 1  shows  a  pointwise  elect  rostatic  approximation  to  the  continuous 
case  using  3  1  point  charges  uniformly  dist  ributed  over  the  b  reast  phantom  surface  as 
indicated  in  Figure  11(b).  The  magnitudes  of  the  point  charges  are  chosen  so  that  they 
equal  the  value  of  the  continuous  charge  dist  ribution  at  that  point .  This  would  serve, 
for  example,  as  a  good  initial  guess  for  a  numerical  optimizat  ion  procedure  which 
would  give  the  optim  um  charge  distributi  on  that  provide  the  m  ost  uniform  electric 
field  magnitude  in  the  breast.  Figure  1 1(c)  shows  the  (unoptimized)  electrostatic  field 
intensity  over  several  planes  in  the  region  to  be  imaged.  It  is  seen  that  even  for  this 
simple  case  that  the  field  intensity  varies  by  less  than  1  dB  in  the  central  region  to  b  e 
imaged.  This  approa  ch  is  cu  rrently  being  deve  loped  an  alytically  as  w  ell  as  b  eing 
modeled  using  two  commercially  available  software  packages. 


Figure  1 1 :  Simulated  uniform  electric  field  (quasi  TEM)  for  a  hemispherical  shell 

excited  by  small  patches. 
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II.4  Experiment 


(a)  Experimental  Setup 

Several  attempts  were  m  ade  to  increase  th  e  electromagnetic  field  in  the  phantom  to  a 
level  sufficient  to  measure  the  resulting  acoustic  signal.  The  ac  oustic  transducer  was 
placed  as  close  to  the  p  hantom  as  possible.  The  specific  acoustic  transducer  used  was 
a  hydrophone,  since  these  are  eh  aracterized  as  having  a  high  sensitivity  over  a  broad 
acoustic  frequency  band.  The  experim  ental  setup  is  illus  trated  in  Fig  ure  12.  The 
exciting  RF  waveguide  aperture  is  extended  to  the  center  of  the  tank  so  a  s  to  increase 
electromagnetic  field  intensity  in  the  vicinity  of  the  phantom. 


Tank 

Filled  with  water 


Figure  12:  Experimental  setup. 

Much  tim  e  has  been  in  vested  into  the  construction  of  appropriate  phantom  s.  The 
current  generation  of  phantom  s  is  bei  ng  constructed  in  c  ooperation  with  the 
Department  of  Biomedical  Engineering  at  the  University  of  Florida.  Earlier  version  of 
the  phantoms  were  made  from  a  mixture  of  a  saline  solution  with  a  gelling  agent,  and 
needed  to  be  held  in  a  latex  container  so  as  to  maintain  the  desired  sh  ape.  Surface 
tension  effects  will  inhibit  the  phantom  from  freely  oscill  ating  to  som  e  degree.  The 
newer  generation  of  phantom  s  is  made  in  a  m  anner  that  provides  a  self-supporting 
structure  and  can  be  placed  directly  in  the  water  bath. 

A  typical  baseband  signal  as  viewed  on  a  spectrum  analyzer  is  shown  in  Figure  13. 
The  anticipated  amplitude  of  the  acoustic  pressure  wave  that  would  be  generated  with 
the  present  experim  ental  setup  is  quite  sm  all,  and  RF  interference  ef  fects,  generated 
by  am  plifier  nonlinearities,  coupled  dire  ctly  though  waveguide  and  tank  onto  the 
metal  housing  of  the  transducer  .  This  is  evidenced  by  the  fact  that  these  sam  e  peaks 
exist  even  in  the  absence  of  the  phantom  inside  the  tank.  Multiple  peaks  are  the  result 
of  the  higher  order  harmonics  of  RF  input  si  gnal.  It  has  been  experim  entally  verified 
that  the  clo  ser  the  tran  sducer  is  placed  to  the  phantom ,  and  hence  to  the  exc  iting 
waveguide,  the  stronger  interference  signal  that  is  picked  up. 
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Figure  13:  Typical  baseband  signal. 

Based  on  these  observations,  a  low  pass  fdter  was  added  after  the  transducer  to  reduce 
the  high  frequency  interference.  The  tank  assembly  was  also  placed  inside  a  Faraday 
shield  to  f  urther  reduce  interference  with  the  measurem  ent  equipment.  The  original 
ultrasonic  transducer  (Y303,  Panametrics  NDT)  was  replaced  with  a  hydrophone  (TC 
4035,  Reson)  which  has  a  broad  bandwidth  (f  rom  10  kHz  to  800  kHz)  and  a  flat 
frequency  response.  The  modified  experimental  setup  is  shown  in  Figure  14. 


Figure  14:  Modified  experiment  setup. 

These  m  odifications  resulted  in  a  subs  tantial  reduction  in  the  aforem  entioned 
interference  signals.  Measurement  of  a  ther  mo-acoustic  signal  with  a  dynam  ic  range 
sufficient  f  or  high  -fidelity  im  aging  w  as  still  p  roblematic.  The  low  R  F  pow  er  lev  el 
used  for  the  RF  excitation  serves  as  the  primary  reason  for  dif  Acuities  encountered 
with  sig  nal  m  easurement.  Pow  er  le  vels  us  ed  in  thes  e  exp  eriments  are  nearly  two 
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orders  of  m  agnitude  lower  than  that  of  other  researchers  [2,  9-1  1].  The  heat 

conduction  into  the  water  bath,  as  previously  discussed,  also  served  to  m  ask  the 
thermo-acoustic  s  ignal.  A  coustic  dissipa  tion  m  echanisms  add  f  urther  dif  Acuities, 
making  it  harder  to  isolate  the  acoustic  pressure  wave  generated  by  the  phantom  from 
other  effects. 


(b)  Experimental  Results  and  Analysis 

The  first  step  in  m  easuring  any  thermoacoustic  resonance  is  the  determination  of  the 
noise  floor  of  the  overall  m  easuring  system,  including  the  equivalent  noise  pressure 
present  at  the  transdu  cer.  Experim  ental  de  termination  of  the  tran  sducer  sens  itivity 
[12],  in  conjunction  with  the  noise  floor  of  the  overall  measurement  system  gives  the 
equivalent  noise  pressure  as  a  function  of  frequency  shown  in  Figure  15.  It  is 
observed  that  the  equivalent  noise  pressure  is  greater  th  an  the  estim  ated  acoustic 
pressure  signal  obtained  from  simulations  .  M  itigation  of  this  problem  requires  a 
lowering  of  the  overall  noise  floor  and/or  increase  in  RF  power  level. 
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Figure  15:  Equivalent  noise  pressure  of  overall  measurement  system. 

It  can  be  sated  with  a  high  degree  of  confidence  that  the  single  most  significant  factor 
resulting  in  the  dif  Acuity  to  d  etect  a  thermo-acoustic  s  ignal  is  in  sufficient  R  F 
excitation  (power).  Other  researchers  have  succeeded  in  obtaining  an  acoustic  signa  1 
used  a  puls  ed  microwave  source  with  peak  power  ranging  from  10  kW  to20kW  [2, 
9-11].  Therefore,  in  order  to  obtain  the  therm  o-acoustic  signal  n  ecessary  for 
generating  im  ages,  a  high  power  pulsed  microwave  source  has  been  obtained. 
Presently  a  pulsed  source  capable  of  s  upplying  0.75  m  icrosecond  pulses  at  a  pulse 
repetition  rate  of  1  100  Hz  and  a  peak  power  up  to  12  OkW  is  cu  rrently  being 

integrated  into  the  ex  perimental  set  up.  (Radio-Resea  rch  Instrum  ent  Co,  P/N 
12- 1-21  MOD).  This  new  experimental  setup  is  shown  in  F  igure  16.  Parts  acquisition 
and  construction  of  this  system  is  currently  underway. 
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Tank  filled  with  water 


Figure  16:  New  experiment  setup. 


With  the  use  of  a  one  m  icrosecond  pulse  d  m  icrowave  source  with  a  peak  power 
ranging  from  20  kW  to  120  kW  ,  it  is  fully  ex  pected  that  the  th  ermo-acoustic  signal 
should  be  found  using  an  oscilloscope.  Th  ere  are  m  any  open  research  areas  where 
improvement  of  the  i  maging  system  performance  can  be  realized,  and  th  e  use  of  the 
high-power  pulsed  source  will  now  allow  investigation  on  these  areas.  Firstly,  the  RF 
aperture  must  be  engineered  to  generate  an  appropriately  conditioned  excitation  signal 
in  the  ph  antom.  Secondly ,  w  ith  th  e  high  pow  er  lev  els  n  ow  availab  le,  nonlin  ear 
phenomenon  can  be  explored.  Since  the  input  power  can  be  high,  the  phantom  itself 
may  exhibit  nonlinear  ef  fects  which  m  ay  be  exploited.  Thirdly,  the  data  acquisition 
system  using  high  speed  ADC’ s  and  a  FPGA  development  platform  can  be  optimized 
for  optimal  signal  integrity  for  imaging  applications.  One  common  solution  would  be 
to  use  of  f-the-shelf  digitizer  cards  to  implement  the  DAQ  sy  stem,  and  t  hese  can  get 
expensive  since  many  channels  are  needed.  However,  one  high  speed  ADC  and  FPGA 
development  board  c  an  be  used  to  im  plement  a  m  ultiplexer-based  d  ata  acqu  isition 
system,  greatly  redu  cing  cost.  Finally  ,  the  aco  ustic  transd  ucers  w  ill  ultim  ately  b  e 
replaced  with  an  array  of  MEMS  sensors  in  an  optical  configuration. 


II.5  Adaptive  Image  Formation  Algorithms 

Developing  accurate  and  robust  im  age  reconstruction  methods  is  o  ne  of  the  key 
challenges  encountered  in  T  AI.  Various  image  reconstruction  algorithms  have  been 
developed  for  T  AI.  By  using  Radon  tran  sformation  on  the  T  AI  data  function, 
reflectivity  tom  ography  reconstruction  al  gorithms  can  be  used  for  T  AI  im  age 
reconstruction  [13].  Exact  in  verse  solutions  have  been  found  for  dif  ferent  scanning 
geometries  in  both  the  frequency  dom  ain  [14,  15]  and  the  tim  e  domain  [16,  17]. 
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Approximate  reconstruction  algorithm  s,  such  as  the  tim  e-domain  Delay-and  -Sum 
(DAS)  beamforming  method  [18,  19]  and  the  optimal  statistical  approach  [20],  have 
also  been  proposed.  However,  a  common  assumption  of  these  existing  methods  is  that 
the  surrounding  tissue  is  acoustically  homogeneous.  This  approximation  is  inadequate 
in  many  medical  imaging  applications.  According  to  previous  studies,  the  sound  speed 
in  hum  an  fern  ale  breast  varies  wide  ly  from  1430m  /s  to  1570m  /s  around  the 
commonly  assum  ed  speed  ofl510m  /s  [  21,  22],  The  heterogeneous  acoustic 
properties  of  biological  tissues  cause  am  plitude  and  phase  distorti  ons  in  the  recorded 
acoustic  signals,  which  can  result  in  significant  degradations  in  imaging  quality. 

Four  robust  and  adaptive  im  age  formation  algorithms,  named  as  adaptive  and  robust 
methods  of  reconstruc  tion  (ARMOR),  multifrequency  adaptive  and  robust  technique 
(MART),  autom  atic  multifrequ  ency  adaptiv  e  and  robust  techniqu  e  (AMART),  and 
iterative  ad  aptive  app  roach  (IAA),  have  been  developed  and  applied  to  the  T  AI 
system. 

ARMOR  is  based  on  robust  Capon  beamfor  ming  (RCB)  [2  3].  This  technique  can  be 
used  to  m  itigate  the  am  plitude  and  phase  disto  rtion  proble  ms  in  T  AI  by  allow  ing 
certain  unc  ertainties.  Sp  ecifically,  in  the  first  step  of  ARMOR,  RCB  is  used  for 
waveform  estim  ation  by  treating  the  am  plitude  distortion  with  a  n  uncer  tainty 
parameter.  In  the  secon  d  step  of  ARMOR,  a  sim  pie  yet  ef  fective  peak  search  ing 
method  is  used  for  phase  distortion  correc  tion.  Com  pared  with  other  ener  gy-  or 
amplitude-based  response  intensity  estimation  methods,  peak  searching  can  be  used  to 
improve  image  quality  with  little  additional  computational  costs.  Moreover,  since  the 
acoustic  pulse  is  usually  bipolar  :  a  positive  pe  ak,  corresponding  to  the  compression 
pulse,  and  a  negative  peak,  corresponding  to  the  rarefact  ion  pulse,  w  e  can  further 
enhance  the  im  age  contras  t  in  T  AI  by  using  the  peak-to-peak  dif  ference  as  the 
response  intensity  for  a  focal  point. 

Instead  of  a  single  f  requency  sour  ce  a  m  ultiple  f  requency  source  is  em  ployed  in 
MART.  MART  can  of  fer  higher  signal-to-ra  tio  (SNR)  and  higher  im  aging  contrast 
than  its  single  frequency  counterpart,  which  we  refer  to  as  the  single-frequency 
adaptive  and  robust  technique  (SART),  since  much  more  information  about  the  human 
breast  can  b  e  harvested  from  the  multiple  frequencies.  Furthermore,  the  inte  rference 
due  to  inhomogeneous  breast  tissue  can  be  suppressed  more  effectively  since  m  ore 
information  about  the  breast  tissue  can  be  used  by  the  RCB  algorithm  .  MART  is  a 
three-stage  tim  e-domain  signal  pro  cessing  a  lgorithm.  In  S  tage  I,  R  CB  is  used  to 
estimate  th  e  therm  al  ac  oustic  resp  onses  f  rom  the  f  ocal  po  ints  w  ithin  the  bre  ast  f  or 
each  stimulating  frequency.  Then  in  Stage  II,  a  scalar  acoustic  waveform  at  each  focal 
point  is  estim  ated  based  on  the  response  estim  ates  for  all  frequencies  from  Stage  I. 
Finally,  in  Stage  III,  the  positive  peak  and  the  negative  peak  of  the  estimated  acoustic 
waveform  at  each  g  rid  location  are  determ  ined,  and  the  p  eak-to-peak  dif  ference  is 
computed  and  referred  to  as  the  image  intensity. 

The  above  two  algorithm  s,  as  well  as  m  ost  of  the  existing  robust  schem  es,  are  us  er 
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parameter  dependent  and  it  may  not  be  a  simple  task  to  determine  the  user  parameters 
in  practice.  Therefore,  user  param  eter-free  rob  ust  adaptive  approaches  ,  including  a 
shrinkage-based  general  linear  combination  (GLC)  algorithm,  are  desirable  [24,  25]. 

We  have  proposed  an  a  utomatic  (i.  e.,  user  param  eter  free)  m  ultifrequency  adaptiv  e 
and  robust  technique  (AMART)  based  on  GLC  for  TAI  to  achieve  high  resolution  and 
good  interference  suppression  capability.  AMART  is  a  three-stage  imaging  algorithm. 
Specifically,  in  the  firs  t  stage  of  AMAR  T,  GLC  is  used  to  estim  ate  the  the  rmal 
acoustic  responses  from  the  grid  points  within  the  breast  for  each  stim  ulating 
frequency.  Based  on  these  estimates,  a  scalar  acoustic  waveform  at  each  grid  point  is 
estimated  via  GLC  at  the  second  stage.  At  th  e  final  stage,  the  energy  of  the  estimated 
acoustic  wavefor  m  at  each  grid  po  int  is  com  puted  and  referred  to  as  the  im  age 
intensity. 
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Figure  17:  Breast  model. 

To  validate  the  ef  fectiveness  of  the  propos  ed  algorithms,  we  have  developed  a  2-D 
inhomogeneous  breast  model,  as  shown  in  Figure  17.  This  breast  model  includes  skin, 
breast  fatty  tissues,  glandular  tissues,  and  the  chest  wall.  Small  tumors  are  set  b  elow 
the  skin.  The  f  inite-difference  time-domain  (FDTD)  method  is  us  ed  to  simulate  the 
electromagnetic  field  inside  the  breast  ti  ssues  [26,  27].  The  specific  absorption  rate 
(SAR)  distribution  is  calculated  based  on  the  simulated  electromagnetic  field  [28,  29]. 
Then  FDTD  is  used  aga  in  to  simulate  the  propagation  of  the  thermal  acoustic  waves 
[30,  3 1].  In  the  followin  g  example,  the  th  ermal  acoustic  sig  nals  are  sim  ulated  based 
on  the  aforem  entioned  2-D  m  odel.  Multip  le  stim  ulating  frequencie  s  from  200-800 
MHz  with  frequency  step  100  MHz  are  used  for  MART.  Two  small  1.5-mm-diameter 
tumors  are  set  insid e  the  breast  model.  Their  locations  are  at  (X=70  mm,  Y=60  mm) 
and  (X=75  mm,  Y=62.5  mm  ).  The  distanc  e  be  tween  the  tw  o  tumors  is  4  m  m.  For 
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comparison  purposes,  the  DAS  m  ethod  is  a  pplied  to  the  sam  e  data  set.  The 
reconstructed  im  ages  are  shown  in  Figur  e  18.  Figure  18(a)  is  the  im  aging  result 
obtained  by  DAS.  The  DAS  image  contains  much  clutter  and  cannot  show  the  tumors 
clearly.  Figures  18(b)  and  18(c)  show  th  e  imaging  results  o  btained  by  A  RMOR  and 
MART.  The  two  turn  ors  are  seen  clearly  in  the  ARMOR  and  MAR  T  images,  and  the 
sizes  and  th  e  locations  of  the  two  turn  ors  are  accurate.  Ho  wever,  the  p  erformance  of 
ARMOR  is  a  little  wors  e  than  that  of  MAR  T  because  so  me  clutter  sh  ow  up  in  the 
ARMOR  images.  Both  ARMOR  and  MART  need  a  user  parameter,  which  is  fixed  at 
0.3N,  where  N  is  the  number  of  the  receiv  er  elements.  Figure  18(d)  is  the  im  aging 
result  obtained  by  the  user  param  eter  free  algorithm,  AMART.  The  image  shows  the 
sizes  and  location  s  of  the  two  turn  ors  accu  rately.  The  advantage  of  AMAR  T, 
compared  with  ARMOR  and  MART,  is  that  it  is  user  parameter  free  and  have  easy  to 
use. 
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Figure  18:  Reconstructed  images,  (a)  DAS,  (b)  ARMOR,  (c)  MART,  and  (d)  AMART. 

Recently,  a  weighted  least  squ  ares-based  non-parametric  and  user  param  eter-free 
iterative  adaptive  approach  (IAA)  [32]  wa  s  proposed  in  array  processing  and  other 
sensing  applications.  IAA  can  work  well  with  few  snapshots  (even  one),  uncorrelated, 
partially  correlated,  and  coherent  so  urces,  and  arbitrary  array  geom  etries.  However, 
due  to  the  properties  of  T  AI,  including  wideband  signals  and  near -field  environment, 
we  cannot  apply  IAA  directly,  since  it  is  de  signed  for  narrowband  signals  originally. 
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We  extend  the  IAA  approach  to  the  wideband  scenario  by  applying  Fourier  transform 
to  the  tim  e-domain  array  output  to  transf  orm  the  wideband  data  into  narrowband 
frequency  bins.  Then  IAA  can  be  applied  to  each  frequenc  y  bin  to  estimate  the  signal 
spectral  distribution  and  hence  estim  ate  the  signal  waveform  and  the  backscattered 
energy.  Figure  19  is  the  im  aging  result  obtained  by  IA  A.  The  breast  m  odel  used  for 
the  simulation  is  the  sam  e  as  the  one  shown  in  Figure  17.  The  early  tim  e  response 
from  the  skin  is  rem  oved.  To  obtain  the  signals,  we  perfor  m  the  simulation  twice  at 
each  stimulating  frequency,  with  and  without  the  tumor,  and  record  the  acoustic  data. 
The  dif  ference  of  the  two  received  signals  is  ref  erred  to  as  the  th  ermal  acoustic 
response  only  from  the  turn  or.  The  t  wo  tumors  are  clearly  shown  in  the  IAA  image, 
and  the  sizes  and  the  locations  of  the  two  tumors  are  accurate. 
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Figure  19:  Reconstructed  image  of  IAA. 


25 


III.  Key  Research  Accomplishments 


•  Phantom  was  developed  and  properties  as  a  function  of  aging  were  explored. 

•  Equipment  needed  for  excitation  and  elect  ric  field  m  easurements  were  identified 
and  purchased  (please  see  Table  1  for  details). 


Vendor 

Model  number 

Equipment  purchased 

Qt. 

Panametrics 

V303-SU  Ultrasoni 

:  Immersion  Transducer 

2 

BCU-74-6W 

Cables.  W  aterproof  &BNC  to  UHF  .6’. 

RG174/U 

2 

BCU-58-10W 

Cable.  W  aterproof  &BNC  to  UHF  .10  ’. 

RG58/U 

1 

5662  Ultrasonic 

Preamplifier 

1 

National 

Instruments 

763000-01  Power 

Cord 

1 

778644-01 

NIPXI-1045  Front  RackMoun  tKit  for  19" 

Rack 

1 

778644-02 

NIPXI-1045  Rear  Rack  Mount  Kit  f  or  19" 

Rack 

1 

778645-01 

NIPXI-1045  18-Slot  3U  PXI  Chassi  s  with 

Universal  AC  Power  Supply 

1 

779505-03 

NI  PXI-PCIe8361,  MXI-Express,  1  Port  PCIe, 

3  m  Cable 

1 

960597-18 

PXI  18-Slot  Factory  I  nstallation  Servi  ce  and 
Extended  Warranty 

1 

778739-01 

NI  PXI-2529  High  Density  Multiconfi  guration 

Matrix 

1 

NIPXI-5122 

Dual  100  MS/s,  14-Bit  Digitizer  w/  anti-alias 
filters  &  8  MB/ch 

16 

778840-01 

NI  TB-2634  Configures  the  NI  PXI-2529  High 
Density  Matrix 

1 

779079-02 

NI  PXI-5671  RF  Signal  Generator 

1 

Shippi 

ng 

1 

TMR 

Plexiglass  Cylinder 

1 

Penn 

Engineering 

1452-4A  W 

aveguide 

1 

6352-5 

Silicon  Gasket  (price  included  in  waveguide) 

Bruce  Carst  en 

Assoc. 

EG01-12"  Probes 

2 

EFP  200R 

Calibrated  radial  probe 

1 

EFP  200A 

Calibrated  axial  probe 

1 

Mini-Circuits  ZI 

IL-100W  -52-S 

Broadband  Amplifier 

1 

Acopian  W530IV 

1T1  3 

Power  supply 

1 

Table  1 .  Equipment  purchased. 
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•  Experimental  setup  to  detect  the  thermal  acoustic  waves  from  the  simulated  tumor 
was  formulated. 

•  A  3-D  electromagnetic  stimulation  system  is  being  simulated  and  developed. 

•  Robust  and  adaptive  im  age  for  mation  algorithms,  including  ARMOR,  MAR  T, 
AMART,  and  IAA,  were  developed. 

•  A  2-D  inhomogeneous  breast  m  odel,  whic  h  inc  ludes  sk  in,  breast  f  atty  tiss  ues, 
glandular  tissues,  and  the  chest  wall,  wa  s  developed  to  validate  the  ef  fectiveness 
of  the  signal  processing  algorithms. 
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VI.  Conclusions 


Numerical  analysis  about  the  electro-acoustic  transduction  of  the  phantom  has  been 
investigated,  based  on  the  Di  ebold’s  research  [5],  in  order  to  g  et  a  more  clear 
understanding  about  the  entire  picture.  The  m  ajor  factors  determining  the  frequency 
response  of  electro-acoustic  transduction  ha  ve  been  pointed  out,  and  som  e  material 
parameters  are  being  measured  based  on  specifically  designed  experiments. 

The  entire  excitation  and  de  tection  system  has  been  set  up,  including  m  odulated 
signal  source,  high  power  amplifier  (100  W),  low  noise  amplifier  and  data  acquisition 
system.  Relative  programming  has  been  com  pleted  and  tested  concerning  the  control 
of  the  entire  system  .  The  noi  se  floor  of  the  m  easurement  system  has  been  obtained 
and  referred  to  the  input  of  transducer  ,  which  provides  the  clear  com  parison  between 
the  signal  level  simulated  and  noise  floor  measured.  Phantoms  are  being  redesigned  to 
gain  m  ore  conductivity  and  hi  gher  resonant  frequency  .  Ot  her  m  odifications  of  the 
system  are  being  conducted  in  order  to  enla  rge  the  signal  level  a  nd  hence,  get  higher 
SNR. 

Four  robust  and  adaptive  im  age  form  ation  algorithm  s,  ARMOR,  MAR  T,  AM  ART, 
and  IAA,  have  been  developed  for  the  T  AI  system.  The  excellent  perform  ance  with 
high  resolution  and  good  inte  rference  suppression  capabilit  y  of  these  algorithm  s  has 
been  demonstrated  based  on  a  2-D  breas  t  model.  Moreover,  the  new  AMAR  T  and 
IAA  algorithms  avoid  the  need  to  s  pecify  any  user  param  eters  and  hence  are  easy  to 
use  in  practice. 
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Abstract — In  this  paper,  we  present  new  adaptive  and  robust 
methods  of  reconstruction  (ARMOR)  for  thermoacoustic  tomog¬ 
raphy  (TAT),  and  study  their  performances  for  breast  cancer  detec¬ 
tion.  TAT  is  an  emerging  medical  imaging  technique  that  combines 
the  merits  of  high  contrast  due  to  electromagnetic  or  laser  stim¬ 
ulation  and  high  resolution  offered  by  thermal  acoustic  imaging. 
The  current  image  reconstruction  methods  used  for  TAT,  such  as 
the  delay-and-sum  (DAS)  approach,  are  data-independent  and  suf¬ 
fer  from  low-resolution,  high  sidelobe  levels,  and  poor  interference 
rejection  capabilities.  The  data-adaptive  ARMOR  can  have  much 
better  resolution  and  much  better  interference  rejection  capabili¬ 
ties  than  their  data-independent  counterparts.  By  allowing  certain 
uncertainties,  ARMOR  can  be  used  to  mitigate  the  amplitude  and 
phase  distortion  problems  encountered  in  TAT.  The  excellent  per¬ 
formance  of  ARMOR  is  demonstrated  using  both  simulated  and 
experimentally  measured  data. 

Index  Terms — Array  signal  processing,  biomedical  acoustic 
imaging,  robustness. 

I.  Introduction 

THERMOACOUSTIC  tomography  (TAT),  the  earliest  in¬ 
vestigation  of  which  dates  back  to  the  1980s  [1],  has  re¬ 
cently  attracted  much  interest  with  its  great  promise  in  a  wide 
span  of  biomedical  applications  (see,  e.g.,  [2]-[4]).  Its  physical 
basis  lies  in  the  contrast  of  the  radiation  absorption  rate  among 
different  biological  tissues.  Due  to  the  thermoacoustic  effect, 
when  a  short  electromagnetic  pulse  (e.g.,  microwave  or  laser) 
is  absorbed  by  the  tissue,  the  heating  results  in  expansion  that 
generates  acoustic  signals.  In  TAT,  an  image  of  the  tissue  ab¬ 
sorption  properties  is  reconstructed  from  the  recorded  thermoa¬ 
coustic  signals.  Such  an  image  may  reveal  the  physiological  and 
pathological  status  of  the  tissue,  which  can  be  useful  in  many 
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applications  including  breast  cancer  detection  [5].  Compared 
with  microwave  imaging  and  ultrasound  imaging,  TAT  com¬ 
bines  their  merits  and  possesses  both  fine  imaging  resolution 
and  good  spatial  contrast  properties  [4] . 

Developing  accurate  and  robust  image  reconstruction  meth¬ 
ods  is  one  of  the  key  challenges  encountered  in  TAT.  Various 
image  reconstruction  algorithms  have  been  developed  for  TAT. 
By  using  Radon  transformation  on  the  TAT  data  function,  re¬ 
flectivity  tomography  reconstruction  algorithms  can  be  used  for 
TAT  image  reconstruction  [6].  Exact  inverse  solutions  have  been 
found  for  different  scanning  geometries  in  both  the  frequency 
domain  [7],  [8]  and  the  time  domain  [9],  [10].  Approximate 
reconstruction  algorithms,  such  as  the  time-domain  delay-and- 
sum  (DAS)  beamforming  method  [11],  [12]  and  the  optimal 
statistical  approach  [13],  have  also  been  proposed.  However, 
a  common  assumption  of  these  existing  methods  is  that  the 
surrounding  tissue  is  acoustically  homogeneous.  This  approx¬ 
imation  is  inadequate  in  many  medical  imaging  applications. 
According  to  previous  studies,  the  sound  speed  in  human  fe¬ 
male  breast  varies  widely  from  1430  to  1570  m/s  around  the 
commonly  assumed  speed  of  15 10  m/s  [14],  [15].  The  heteroge¬ 
neous  acoustic  properties  of  biological  tissues  cause  amplitude 
and  phase  distortions  in  the  recorded  acoustic  signals,  which 
can  result  in  significant  degradation  in  imaging  quality  [16]. 

In  ultrasound  tomography  (UT),  wavefront  distortion  due  to 
heterogeneity  of  biological  tissue  has  been  studied  extensively. 
Various  wavefront  correction  methods  have  been  proposed  [17]. 
However,  they  are  not  highly  effective  at  correcting  severe  am¬ 
plitude  distortions  [18],  and  they  usually  involve  complicated 
procedures.  The  problem  in  TAT  is  somewhat  different  from  that 
in  UT.  In  the  breast  UT,  the  amplitude  distortion  caused  by  re¬ 
fraction  is  more  problematic  than  the  phase  distortion  induced  by 
acoustic  speed  variation.  In  TAT,  however,  even  for  the  biologi¬ 
cal  tissue,  such  as  the  breast  tissue,  with  a  relatively  weak  hetero¬ 
geneity,  phase  distortion  dominates  amplitude  distortion  [16]. 
These  unique  features  suggest  that  new  adaptive  and  robust 
imaging  techniques  should  be  designed  especially  for  TAT. 

Time-domain  approximate  reconstruction  algorithms,  such 
as  the  DAS  (weighted  or  unweighted)  type  of  data-independent 
approaches  have  various  applications  in  medical  imaging.  They 
need  little  prior  information  on  the  tissue  for  image  reconstruc¬ 
tion  and  can  be  fast  and  simple  to  implement  to  process  the 
wideband  acoustic  signals.  Although  not  based  on  the  exact 
solution,  they  provide  similar  image  qualities  to  those  of  the  ex¬ 
act  reconstruction  algorithms.  However,  these  data-independent 
methods  tend  to  suffer  from  poor  resolution  and  high-sidelobe- 
level  problems.  Data-adaptive  approaches,  such  as  the  recently 
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introduced  robust  Capon  beamforming  (RCB)  method  [19],  can 
have  much  better  resolution  and  much  better  interference  rejec¬ 
tion  capability  than  their  data-independent  counterparts. 

We  propose  adaptive  and  robust  methods  of  reconstruction 
(ARMOR)  based  on  RCB  for  TAT.  ARMOR  can  be  used  to 
mitigate  the  amplitude  and  phase  distortion  problems  in  TAT 
by  allowing  certain  uncertainties.  Specifically,  in  the  first  step 
of  ARMOR,  RCB  is  used  for  waveform  estimation  by  treating 
the  amplitude  distortion  with  an  uncertainty  parameter.  In  the 
second  step  of  ARMOR,  a  simple,  yet  effective,  peak  searching 
method  is  used  for  phase  distortion  correction.  Compared  with 
other  energy-  or  amplitude-based  response  intensity  estimation 
methods,  peak  searching  can  be  used  to  improve  image  quality 
with  little  additional  computational  costs.  Moreover,  since  the 
acoustic  pulse  is  usually  bipolar:  a  positive  peak,  corresponding 
to  the  compression  pulse,  and  a  negative  peak,  corresponding 
to  the  rarefaction  pulse  [1 1],  we  can  further  enhance  the  image 
contrast  in  TAT  by  using  the  peak-to-peak  difference  as  the 
response  intensity  for  a  focal  point.  We  will  demonstrate  the 
excellent  performance  of  ARMOR  by  using  both  data  simulated 
on  a  2-D  breast  model  and  data  experimentally  measured  from 
mastectomy  specimens. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  II,  we  formulate  the  problem  of  interest.  Sections  III-V 
describe  the  first,  second,  and  third  steps  of  ARMOR,  respec¬ 
tively.  Examples  based  on  simulated  and  real-world  experimen¬ 
tal  data  are  presented  in  Section  VI.  Finally,  Section  VII  provides 
the  conclusions. 


II.  Problem  Formulation 

Consider  a  TAT  imaging  system,  as  shown  in  Fig.  1(a).  A 
stimulating  electromagnetic  (laser  or  microwave)  pulse  is  ab¬ 
sorbed  by  the  biological  tissue  under  testing,  which  causes  a 
sudden  heat  change  (of  the  order  of  10-4  °C  [20]).  Due  to  the 
thermoacoustic  effect,  an  acoustic  pulse  is  generated  that  can 
be  recorded  by  an  ultrasonic  transducer  array.  The  transducer 
array  may  be  a  real  aperture  array  or  a  synthetic  aperture  array 
formed  by  rotating  a  sensor  around  the  tissue  and  recording  the 
acoustic  waves  at  different  locations.  We  assume  that  the  num¬ 
ber  of  transducers  in  the  array  (or  in  the  synthetic  aperture  array 
case,  the  number  of  transducer  data  acquisition  locations)  is  M. 
Each  transducer  is  assumed  to  be  omnidirectional;  mutual  cou¬ 
plings  among  the  transducers  are  not  considered  in  our  model 
as  they  can  be  tolerated  by  our  robust  algorithms  to  a  certain 
extent.  The  recorded  acoustic  signals  are  sufficiently  sampled 
and  digitized  and  a  typical  recorded  pulse  is  shown  in  Fig.  1(b) 
(based  on  the  data  measured  on  the  breast  specimen  II  described 
in  Section  VI). 

The  data  model  for  the  sampled  and  digitized  acoustic  signal 
recorded  by  the  mth  transducer  is  given  by: 

xm(n)  =  sm(n)  +  em(n),  m  -•  1 . M.  (1) 

where  n  is  the  discrete  time  index,  starting  from  to  after  the  ex¬ 
citation  pulse.  The  scalar  sm  ( n )  denotes  the  signal  component, 
which  corresponds  to  the  acoustic  pulse  generated  at  a  focal 
point,  and  em  (n)  is  the  residual  term,  which  includes  unmod¬ 


Fig.  1.  (a)  A  schematic  of  a  2-D  synthetic-aperture-based  TAT  scanning  sys¬ 

tem.  (b)  A  typical  acoustic  pulse  recorded  by  a  transducer  (for  data  measured 
from  breast  specimen  II). 


eled  noise  and  interference  (caused  by  other  sources  within  the 
tissue). 

The  goal  of  ARMOR  is  to  reconstruct  an  image  of  thermoa¬ 
coustic  response  intensity  /( r),  which  is  directly  related  to  the 
absorption  property  of  the  tissue,  from  the  recorded  data  set 
{xm  (n)}.  Herein,  the  (2-D  or  3-D)  vector  r  denotes  the  focal 
point  location  coordinate.  To  form  an  image,  we  scan  the  focal 
point  location  r  to  cover  the  entire  cross  section  of  the  tissue 
(the  transducers  can  acquire  signals  at  different  heights;  for  each 
height,  a  2-D  cross-sectional  image  can  be  reconstructed  and  a 
3-D  image  can  be  formed  from  the  2-D  images).  We  allow  cer¬ 
tain  uncertainties  in  ARMOR  to  deal  with  amplitude  and  phase 
distortions  caused  by  the  background  heterogeneity. 

The  discrete  arrival  time  of  the  pulse  (for  the  mth  transducer) 
can  be  determined  approximately  as 


im(r) 


to_  ||  r  —  rm  1| 
At  Atvo 


(2) 


We  will  omit  the  dependence  of  the  arrival  time  tm  (r)  on  r 
hereafter  for  notational  simplicity.  Here,  At  is  the  sampling 
interval,  and  the  3-D  vector  rm  denotes  the  location  of  the  mth 
transducer.  The  sound  speed  Vo  is  chosen  to  be  the  average 
sound  speed  of  the  biological  tissue  under  interrogation.  The 
notation  ||x||  denotes  the  Euclidean  norm  of  x,  and  [y\  stands 
for  rounding  to  the  greatest  integer  less  than  y.  The  second  term 
in  (2)  represents  the  time-of- flight  between  the  focal  point  and 
the  mth  transducer. 

The  signal  components  {sm{p)}m=i  are  approximately 
scaled  and  shifted  versions  of  a  nominal  waveform  s(t)  at  the 
source 


exp  (— a\\r 


where  a  is  the  attenuation  coefficient  in  Nepers/m.  In  TAT, 
the  major  frequency  components  of  the  acoustic  signals  take  a 
relatively  narrow  band,  and  are  usually  lower  than  those  in  UT 
[16].  Hence,  we  can  approximate  a  as  a  frequency-independent 
constant. 

We  preprocess  the  data  to  time  delay  all  the  signals  from  the 
focal  point  r  and  compensate  for  the  loss  in  amplitude  due  to 
propagation  decay.  Fet  ym  (n)  denote  the  signal  after  prepro¬ 
cessing  to  backpropagate  the  detected  signal  to  the  source 


ym{n)  =  exp(a||r  -  rm||)  •  ||r-rm||  • xm(n  +  tm ).  (4) 
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Then,  the  received  vector  data  model  can  be  written  as 

y(n)  =  a0s(n)  +  e(n),  n=  (5) 

where  a0  is  the  corresponding  steering  vector,  which  is  approxi- 
mately  equal  to  a  =  [1, . . . ,  l]r,  y(n)  =  [yi(n), . . .  ,yM  (n)]T, 
e(n)  represents  the  noise  and  interference  term  after  preprocess¬ 
ing,  and  (*)T  denotes  the  transpose.  Here,  we  define  the  time 
interval  of  interests  for  the  signal  y(t)  to  be  from  —TV  to  TV, 
which  means  that  we  only  take  N  samples  before  and  after  the 
approximate  arrival  time  given  in  (2)  for  the  focal  point  at  r.  The 
value  of  N  should  be  chosen  large  enough  so  that  the  interval 
from  —N  to  N  covers  the  expected  signal  duration  in  the  region 
of  interest. 

In  reality,  both  the  amplitude  and  the  phase  (or  pulse  arrival 
time)  of  the  acoustic  pulse  will  be  distorted.  A  major  cause  for 
these  distortions  is  the  acoustically  heterogeneous  background. 
Amplitude  distortion  is  mainly  due  to  the  interferences  caused 
by  multipath,  which  is  inevitable  in  the  heterogeneous  medium: 
refraction  occurs  due  to  acoustic  speed  mismatch  across  the  tis¬ 
sue  interface;  consequently,  acoustic  pulses  arrived  at  the  trans¬ 
ducer  will  be  via  different  routes  and  interfere  with  each  other. 
On  the  other  hand,  phase  distortion  is  mainly  caused  by  the 
nonuniform  sound  speed.  For  example,  in  human  female  breast, 
the  sound  speed  can  vary  from  1430  to  1570  m/s;  therefore, 
the  actual  arrival  time  will  fluctuate  around  the  approximately 
calculated  time  given  in  (2).  Moreover,  an  inaccurate  estimate 
of  to  (to  is  aligned  with  the  focal  point’s  signal  arrival  time) 
and  the  transducer  calibration  error  may  also  contribute  to  the 
phase  distortion.  Amplitude  and  phase  distortion  will  blur  the 
image,  raise  the  image  background  noise  level,  lower  the  values 
of  the  object  of  interest,  and,  consequently,  decrease  the  image 
contrast  [16]. 

We  mitigate  the  effects  of  these  distortions  by  allowing  ao  to 
belong  to  an  uncertainty  set  centered  at  a  and  by  considering 
the  signal  arriving  within  the  interval  from  —NtoN. 

III.  Step  I  of  ARMOR:  Waveform  Estimation 

The  first  step  of  ARMOR  is  to  estimate  the  waveform  of  the 
acoustic  pulse  generated  by  the  focal  point  at  location  r,  based 
on  the  data  model  in  (5).  It  will  appear  that  we  have  neglected  the 
presence  of  phase  distortion  by  using  this  data  model  in  the  first 
step.  However,  by  allowing  ao  to  be  uncertain,  we  can  tolerate 
some  phase  distortions  as  well.  This  approximation  causes  little 
performance  degradation  to  our  robust  algorithm. 

Covariance-fitting-based  RCB  [21]  is  used  to  first  estimate 
the  steering  vector  ao ,  and  use  the  estimated  ao  to  obtain  an  op¬ 
timal  beamformer  weight  vector  for  pulse  waveform  estimation. 
By  assuming  that  the  true  steering  vector  lies  in  the  vicinity  of 
the  nominal  steering  vector  a,  we  consider  the  following  opti¬ 
mization  problem  [19] 

max  cr2  subject  to  R  — cr2aoao  ^0, 

cr2,a0 

||ao  -  a||2  <  e,  (6) 


where  A  y  0  means  that  the  matrix  A  is  positive  semidefinite, 
cr2  is  the  power  of  the  signal  of  interest,  and 

1  N 

ft  =  2jv  + 1  y(n)yT(n)  (?) 

is  the  sample  covariance  matrix.  The  second  constraint  in  (6)  is 
a  spherical  uncertainty  set;  an  elliptical  uncertainty  set  can  be 
used  instead,  if  a  tighter  constraint  is  desirable  [21]. 

The  parameter  5  in  (6)  determines  the  size  of  the  uncertainty 
set  and  is  a  user  parameter.  To  avoid  the  trivial  solution  of 
ao  =  0,  we  require  that 

£<||a||2.  (8) 


It  can  be  verified  that  the  smaller  the  5,  the  higher  the  resolution 
and  the  stronger  the  ability  of  RCB  to  suppress  an  interference 
that  is  close  to  the  signal  of  interest,  and  that  the  larger  the  5, 
the  more  robust  RCB  will  be  to  tolerate  distortions  and  small- 
sample- size  problems  caused  by  calculating  R  in  (7)  from  a 
finite  number  of  data  vectors  or  snapshots.  When  e  is  close  to 
M,  RCB  will  perform  like  DAS.  To  attain  high  resolution  and 
to  effectively  suppress  interference,  5  should  be  made  as  small 
as  possible.  On  the  other  hand,  the  smaller  the  sample  size  N 
or  the  larger  the  distortions,  the  larger  should  e  be  chosen  [19]. 
Since  the  performance  of  RCB  does  not  depend  very  critically 
on  the  choice  of  6  (as  long  as  it  is  set  to  be  a  “reasonable 
value”)  [21],  such  qualitative  guidelines  are  usually  sufficient 
for  making  a  choice  of  6.  We  will  investigate  the  effect  of  e  in 
Section  VI.  In  our  examples  in  Section  VI,  we  choose  certain 
reasonable  initial  values  for  6,  and  then  make  some  adjustments 
empirically  based  on  image  quality:  making  it  smaller  when  the 
resulting  images  have  low  resolution,  or  making  it  larger  when 
the  image  is  distorted  by  interferences. 

By  using  the  Lagrange  multiplier  method,  the  solution  to  (6) 
is  given  by  [19] 

a0  =  a  —  [I  +  /uR]-1a  (9) 

where  I  is  the  identity  matrix,  and  fi  >  0  is  the  correspond¬ 
ing  Lagrange  multiplier  that  can  be  solved  from  the  following 
equation 

||(I  +  /zR)-1a||2  =  e.  (10) 

Consider  the  eigendecomposition  on  the  sample  covariance  ma¬ 
trix  R 


R  =  urur  (ii) 

where  the  columns  of  U  are  the  eigenvectors  of  R  and  the 
diagonal  matrix  T  consists  of  the  corresponding  eigenvalues 
7i  >  72  >  •  •  •  >  7 m  •  Let  b  =  UT  a,  where  bm  denotes  its  rath 
element.  Then,  (10)  can  be  rewritten  as 


M 

Am)  = 

m  =  1 


\K\ 2 

(1  +  M7m)2 


=  5. 


(12) 


Note  that  C(fjb)  is  a  monotonically  decreasing  function  of  with 
£(0)  >  e  by  (8)  and  lim^oo  C{/i)  =  0  <  5,  which  means  that 
/i  can  be  solved  efficiently,  say,  by  using  the  Newton’s  method 
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(see  [19]  for  more  details).  After  obtaining  the  value  of  yU,  the 
estimate  ao  of  the  actual  steering  vector  ao  is  determined  by  (9). 

Observe  that  there  is  a  “scaling  ambiguity”  in  (6)  by  treating 
both  the  signal  power  a2  and  the  steering  vector  a0  as  un¬ 
knowns  (see  [19]  and  [21]).  The  ambiguity  exists  in  the  sense 
that  (cr2,  ao)  and  (cr2/c,  c1/2 ao)  (for  any  constant  c  >  0)  yield 
the  same  term  cr2  ao  aj .  To  eliminate  this  ambiguity,  we  scale 
the  solution  ao  to  make  its  norm  satisfy  the  following  condition 

||ao||  2=M.  (13) 


(NotethatM=  ||a||2.) 

To  obtain  an  estimate  for  the  signal  waveform  s(n),  we  apply 
a  weight  vector  to  the  preprocessed  signals  {y{n)}^=_N .  The 
weight  vector  is  determined  by  using  the  estimated  steering 
vector  ao  in  the  weight  vector  expression  of  the  standard  Capon 
beamformer  (see,  e.g.,  [19]  and  [21]) 


w  - 

WRCB  -  77172 


R+  il 

ll 

-l 

a0 

a0 

R  +  -I 

L  ^  J 

-1  * 

R 

“  i  A 

-i 

a0 

The  bipolar  acoustic  pulse  has  one  peak  positive  and  another 
negative.  We  determine  the  positive  and  negative  peak  values  as 
follows: 


P+  =  max  < 

[  max  s(n)ol  , 

(18) 

{ne[~, A, A]  J 

P  =  min  < 

min  s(n)  0  >  , 

(19) 

l 

^ne[-A,A]  J 

where  the  searching  range 

[-A,  A]  e  [-N,N] 

is  around  the 

calculated  arrival  time  given 

by  (2).  Here  A  is  a  user  parameter. 

Since  the  peak  searching  is  independent  of  the  particular  wave¬ 
form  estimation  methods,  we  use  s(n)  to  denote  the  waveform 
estimated  by  either  DAS  or  ARMOR. 

The  search  range  is  determined  by  the  difference  between  the 
true  arrival  time  tm  and  the  calculated  arrival  time  tm ,  based 
on  (2).  This  arrival  time  difference  has  been  analyzed  for  breast 
tissue  by  taking  into  account  its  relatively  weak  heterogeneity 
acoustic  property  [16].  An  expression  for  this  difference  is  given 
in  [16]  by 


Note  that  (14)  has  a  diagonal  loading  form,  which  allows  the 
sample  covariance  matrix  to  be  rank-deficient.  The  beamformer 
output  can  be  written  as 

srcbO)  =  WRCBy(n),  n  =  -N, . . . ,  N  (15) 

which  is  the  waveform  estimate  for  the  acoustic  pulse  generated 
at  the  focal  point  at  location  r. 

RCB  can  provide  a  much  better  waveform  estimate  than  the 
conventional  DAS  but  at  a  higher  computational  cost.  For  a 
single  focal  point,  RCB  requires  0(M3)  flops,  which  mainly 
come  from  the  eigendecomposition  of  the  sample  covariance 
matrix  R  [19];  DAS  needs  only  O(M)  flops.  DAS  can  be  used 
as  a  fast  image  reconstruction  method  to  provide  initial  imaging 
results. 

The  weight  vector  used  by  DAS  for  waveform  estimation  is 
wdas  =  a  (16) 

and  the  estimated  waveform  is  given  by 

M 

SDAS (n)  =  WpASy(n)  =  ^  (n),  n  =  -N, . . .  ,N. 

m  =  1 

(17) 


IV.  Step  II  of  ARMOR:  Peak  Searching 

Based  on  the  estimated  waveform  obtained  in  Step  I  for  the 
focal  point  at  location  r,  in  Step  II  of  ARMOR,  we  will  search 
for  the  two  peaks  of  the  bipolar  acoustic  pulse  generated  by 
the  focal  point.  In  a  homogeneous  background,  where  phase 
distortion  is  absent,  we  can  accurately  calculate  the  arrival  time 
of  the  acoustic  pulse  generated  by  the  focal  point  at  location  r  by 
using  (2).  However,  this  is  never  true  in  heterogeneous  biological 
tissues.  It  was  reported  in  [16]  that  when  the  heterogeneity  is 
weak,  such  as  in  the  breast  tissue,  amplitude  distortion  caused 
by  multipath  is  not  severe.  We  can  assume  that  the  original 
peak  remains  a  peak  in  the  waveform  estimated  from  Step  I  of 
ARMOR. 


(r  )  —  tm  tm  OC 


Mr')  -  u0] 


(20) 


where  r'  is  a  point  within  the  line  connecting  the  focal  point 
at  location  r  and  the  rath  transducer  at  location  rm,  and  v(rf) 
is  the  local  sound  speed.  The  higher  order  terms  of  [v(r')  — 
vq\/vq  in  (20)  have  been  ignored.  It  is  reasonable  to  assume 
that  v(rf)  is  Gaus sian-distributed  with  mean  vq  and  variance 
cr2 .  Consequently,  the  arrival  time  difference  is  also  Gaussian- 
distributed  with  zero-mean  and  variance  cr2  oc  cr2  jv\ .  If  we 
choose  A  =  and  the  duration  of  the  acoustic  pulse  is  r,  we 
can  find  the  two  peaks  of  the  pulse  within  the  interval  (— + 
r)  on  the  recorded  signals  with  a  high  probability  of  0.6826. 
This  analysis  is  consistent  with  the  experimental  measurements 
in  [22].  From  our  examples,  we  found  that  a  symmetric  range 
[—A,  A]  around  the  estimated  arrival  time  performs  similarly 
to  the  asymmetric  range  [— A,  A  +  r],  and  we  use  the  former 
since  it  is  easy  to  handle  in  practice.  Also,  we  can  use  similar 
techniques  as  those  in  [22]  to  estimate  <j§  to  find  a  good  searching 
range  for  Step  II  of  ARMOR,  and  to  estimate  r  for  the  energy- 
type  methods,  as  shown  in  our  examples  later. 

There  is  a  tradeoff  in  choosing  the  searching  range.  The  larger 
the  searching  range,  the  higher  the  probability  we  can  find  the 
peaks  of  the  acoustic  pulse  within  the  range.  However,  if  the 
range  is  chosen  too  large,  the  interferences  may  cause  false 
peaks,  and  as  a  consequence,  we  are  more  likely  to  find  a  false 
peak.  In  our  examples  in  Section  VI,  we  choose  the  best  search¬ 
ing  range  empirically  based  on  the  estimated  variance  of  the 
arrival  time  difference  &s . 


V.  Step  III  of  ARMOR:  Intensity  Calculation 

After  estimating  the  waveform  generated  by  the  focal  point  at 
location  r,  we  need  to  obtain  the  response  intensity  based  on  the 
estimated  waveform.  For  the  same  estimated  waveform,  differ¬ 
ent  approaches  can  be  used  to  evaluate  the  focal  point  response 
intensity.  These  approaches  extract  different  information  from 
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the  estimated  waveform  as  the  response  intensity,  and  may  be 
useful  to  physicians  in  different  ways. 

There  are  two  major  types  of  response  intensity  measurement 
approaches:  amplitude-based  and  energy-based.  The  waveform 
peak  values  obtained  in  Step  II  of  ARMOR  can  be  used  for  both 
approaches. 

Conventional  DAS  uses  the  amplitude-based  measure  for  TAT 
imaging  [11],  [12],  with  the  corresponding  response  intensity 
given  by  s(0),  or  equivalently 


M 

ic=s(Q)  =  YJym{Q)  (2D 

m  —  1 


where  the  subscript  “c”  stands  for  “Conventional.” 

The  energy -based  measure,  such  as  the  one  used  in  [23], 
calculates  the  response  intensity  as  follows 


IeI  =  £2(0) 


Vm  (0) 

_m  =  1 


(22) 


where  the  subscript  “ei”  means  “Energy-type  1.” 

The  entire  pulse  energy  has  also  been  used  as  an  intensity 
measure,  such  as  in  the  monostatic  and  multistatic  microwave 
imaging  for  breast  cancer  detection  [24],  [25],  and  the  intensity 
is  given  by 


lE2  =  ±s2(n)  =  Y 


n  —  0 


n  =  0 


M 


~i  2 


Y ym 


,m  =  1 


(23) 


where  the  subscript  “e2”  stands  for  “Energy-type  2.” 

We  can  consider  using  the  peak  value  as  the  response  intensity 
measure  due  to  the  bipolar  nature  of  the  response  at  the  focal 
point 


fP+,  if  |P+|  >  |P-| 

V  =  ’  '  'A  (24) 

L  P  3  otherwise 


where  the  subscript  “p”  stands  for  “Peak,”  with  P+  and  P~ 
defined  in  (18)  and  (19),  respectively.  Herein,  we  keep  the  sign 
of  the  maximum  amplitude  since  the  sign  of  the  peak  may  also 
contain  some  information  about  the  focal  point. 

Peak  searching  maximizes  the  output  signal-to-noise  ratio. 
An  intuitive  explanation  is  that,  given  the  fact  that  the  acoustic 
pulse  is  bipolar  [1 1],  if  we  assume  that  the  residual  term  e(t)  is 
stationary,  or  its  power  is  uniform  over  time,  then  the  signal-to- 
noise  ratio  (SNR)  is  maximized  at  the  (positive  or  negative)  peak 
of  the  acoustic  pulse.  As  a  comparison,  the  conventional  DAS 
(21)  fixes  the  samples  to  be  summed  up  at  the  calculated  arrival 
time.  Due  to  phase  distortions,  the  waveform  at  the  calculated 
time  may  be  far  from  the  peak  value. 

We  can  also  employ  peak-to-peak  difference  as  the  response 
intensity  for  the  focal  point  at  location  r 

Ipp  =  P+  -  P~  >  0  (25) 


where  the  subscript  “pp  ”  denotes  the  “peak-to-peak  difference.” 
Peak-to-peak  difference  has  higher  imaging  contrast  than  peak 
value  measure:  the  peak-to-peak  difference  of  the  bipolar  pulse 
is  approximately  twice  the  absolute  peak  value,  which  means 
that  the  output  signal  power  of  the  former  is  four  times  of  the 


Fig.  2.  2-D  breast  model  in  an  x-y  coordinate  system,  with  a  2-mm-diameter 

tumor  present,  (a)  Model  for  electromagnetic  simulation,  (b)  Model  for  acoustic 
simulation. 


latter;  yet,  the  noise  power  of  the  former  may  be  only  twice 
that  of  the  latter.  Therefore,  the  output  SNR  may  be  doubled 
by  using  the  peak-to-peak  difference  rather  than  the  peak  value. 
Both  peak- value  and  peak-to-peak  difference  measures  belong 
to  the  amplitude-based  measures. 

VI.  Numerical  and  Experimental  Examples 

We  demonstrate  the  performance  of  ARMOR  using  both  nu¬ 
merically  simulated  and  experimentally  measured  TAT  data. 
The  ARMOR  images  are  compared  with  the  DAS  images. 


A.  Numerical  Examples 

We  consider  a  2-D  breast  model,  as  shown  in  Fig.  2.  The  2-D 
breast  model  includes  2-mm  thick  skin,  chest  wall,  as  well  as 
randomly  distributed  fatty  breast  tissues  and  glandular  tissues. 
The  cross  section  of  the  breast  model  is  a  half-circle  with  a  10  cm 
diameter.  In  the  first  numerical  example,  a  2-mm-diameter  tumor 
is  located  at  2.2  cm  below  the  skin  (at  x  =  7.0  cm,  y  =  6.0  cm). 
Fig.  2  shows  the  shape,  dielectric  properties,  and  sound  speed 
variations  of  the  breast  model,  as  well  as  the  tumor  size  and 
location  for  the  first  example.  In  the  second  numerical  example, 
one  large  tumor  (1  cm  in  diameter)  is  located  at  x  =  12  cm, 
y  =  15  cm.  Other  properties  of  the  breast  model  for  the  second 
example  are  the  same  as  those  for  the  first  example. 

To  reduce  the  reflections  from  the  skin,  the  breast  model 
is  immersed  in  a  lossless  liquid  with  permittivity  similar  to 
that  of  the  breast  fatty  tissue.  Seventeen  transducers  (assumed 
omnidirectional)  are  located  on  a  half-circle  10  mm  away  from 
the  skin,  with  uniform  spacing,  to  form  a  real  aperture  array. 

The  dielectric  properties  of  the  breast  tissues  are  assumed  to 
be  Gaussian  random  variables  with  variations  of  ±10%  around 
their  nominal  values.  This  variation  represents  the  upper  bound 
reported  in  the  literature.  The  nominal  values  are  chosen  to  be 
typical  of  those  reported  in  the  literature  [5],  [26],  which  is 
given  in  Table  I  [24] .  The  dielectric  constants  of  glandular  tis¬ 
sues  are  between  er  =  11  and  er  =  15.  The  dispersive  properties 
of  the  fatty  breast  tissue  and  those  of  the  tumor  are  also  consid¬ 
ered  in  the  model.  The  randomly  distributed  breast  fatty  tissues 
and  glandular  tissues  with  variable  dielectric  properties  are  rep¬ 
resentative  of  the  nonhomogeneity  of  the  breast  of  an  actual 
patient. 


Authorized  licensed  use  limited  to:  University  of  Florida.  Downloaded  on  January  29,  2010  at  14:45  from  IEEE  Xplore.  Restrictions  apply. 


2746 


IEEE  TRANSACTIONS  ON  BIOMEDICAL  ENGINEERING,  VOL.  55,  NO.  12,  DECEMBER  2008 


TABLE  I 
Acronyms 


ART 

Adaptive  and  Robust  Methods  Of  Reconstruction 

DAS 

Delay-And-Sum 

FDTD 

Finite  Difference  Time  Domain 

PML 

Perfectly  matched  layer 

RCB 

Robust  Capon  Beamforming 

SNR 

Signal-to-Noise  Ratio 

SAR 

Specific  Absorption 

TAT 

Thermoacoustic  Tompgraphy 

UT 

Ultra-sound  Tomography 

C 

Conventional 

El 

Energy-type  1 

E2 

Energy-type  2 

P 

Peak 

PP 

Peak-to-Peak  difference 

Following  the  report  that  the  breast  tissues  have  a  weak 
acoustic  heterogeneity  [16],  we  model  the  sound  speed  within 
the  breast  as  a  Gaussian  random  variable  with  variation  ±5% 
around  the  assumed  average  sound  speed  of  1500  m/s.  Since 
the  attenuation  coefficient  a  in  (3)  is  small  for  breast  tissue 
(0.75  dB/(MHz-cm))  [15]  and  the  acoustic  signals  are  below 
2  MHz,  we  neglect  the  exponential  attenuation  in  acoustic  wave 
propagation.  Also,  since  the  acoustic  pressure  field  generated 
by  the  thermoacoustic  effect  is  usually  small  [20],  we  do  not 
consider  the  nonlinear  acoustic  effects.  The  probing  microwave 
pulse  used  here  is  a  modulated  rectangular  pulse  with  a  mod¬ 
ulating  frequency  of  800  MHz.  The  duration  of  the  pulse  is 
1  fi s.  More  details  about  the  thermal  acoustic  simulations  are 
given  in  the  Appendix.  In  the  following,  all  the  images  are  dis¬ 
played  on  a  linear  scale,  and  we  will  name  the  imaging  methods 
by  their  waveform  estimation  method  followed  by  the  intensity 
calculation  approach,  such  as  “DAS-C.” 

Note  that  the  skin  also  absorbs  microwave  energy  and  gen¬ 
erates  acoustic  signals.  The  skin  response  is  much  stronger 
than  that  of  the  tumor,  since  the  skin  has  a  much  larger  area 
than  the  tumor  and  the  skin  is  closer  to  the  acoustic  sensors. 
So,  before  applying  the  aforementioned  preprocessing  steps 
and  ARMOR,  we  remove  the  strong  skin  response  using  tech¬ 
niques  similar  to  those  in  [24] .  A  calibration  signal  is  obtained 
as  the  average  of  the  recorded  signals  containing  similar  skin 
response.  Then,  the  calibration  signal  is  subtracted  out  from 
all  recorded  signals  to  remove  the  skin  response  as  much  as 
possible. 

The  searching  range  is  chosen  by  the  guidelines  presented  in 
Section  IV.  To  obtain  a  general  profile  of  the  arrival  time  dif¬ 
ference  caused  by  the  phase  distortion,  we  use  a  simple  method 
similar  to  the  one  used  in  [27].  First,  the  cross-correlation  func¬ 
tions  for  all  the  signals  recorded  by  the  two  adjacent  transducers 
are  obtained.  The  peak  value  of  the  cross-correlation  function 
is  used  to  estimate  the  arrival  time  delay  between  the  signals 
recorded  by  the  adjacent  transducers.  Second,  these  arrival  time 
delays  are  fitted  using  a  fourth-order  polynomial  curve,  which 
is  dominated  by  the  arrival  time  delays  due  to  the  path  length 
differences  in  the  absent  of  phase  distortions.  The  fourth-order 
polynomial  is  used  since  the  delay  caused  by  the  path  length 
difference  should  vary  smoothly  [27].  Fig.  3(a)  shows  the  esti¬ 
mated  arrival  time  delay  and  the  delay  based  on  curve  fitting. 


Fig.  3.  (a)  Comparison  between  the  estimated  and  fitted  arrival  time  delays, 

for  the  simulated  breast  model  with  one  tumor  (the  curves  for  the  two-tumor 
case  are  similar).  Histograms  of  delay  differences,  (b)  Simulated  breast  model 
with  one  tumor,  (c)  Breast  specimen  I.  (d)  Breast  specimen  II. 


Third,  the  delay  difference  between  the  estimated  arrival  time 
delay  and  the  fitted  delay,  or  the  fitting  error,  is  treated  as  the 
arrival  time  distortion  for  the  transducers.  The  standard  devia¬ 
tion  of  the  delay  difference  is  used  to  estimate  as .  Although  the 
accuracy  of  the  cross-correlation  method  is  limited  due  to  false 
peaks  and  jitter  problems,  it  is  sufficient  to  obtain  a  qualitative 
profile  for  as . 

Fig.  3  gives  the  histogram  of  the  delay  difference  for  all  the 
cases  that  we  considered  herein.  For  the  simulated  example, 
the  standard  deviation  of  the  delay  difference  is  4.5,  which 
indicates  a  weak  phase  distortion  in  the  breast  model.  We  set 
an  initial  value  for  A,  based  on  the  estimated  as,  and  adjust 
the  length  of  the  searching  range  to  achieve  the  best  imaging 
result. 

To  estimate  the  pulse  duration  f  (used  in  DAS-E2  and 
RCB-E2),  we  select  several  typical  signals  (with  clear  peaks) 
and  take  the  average  of  their  pulse  durations.  In  practice,  the 
acoustic  pulse  duration  is  determined  by  the  probing  pulse  du¬ 
ration,  size  and  shape  of  the  tumor,  as  well  as  the  transducer 
response. 

Fig.  4  shows  the  images  for  the  simulated  breast  model  with 
one  2-mm  diameter  tumor  formed  using  ARMOR  and  DAS. 
The  tumor  response  is  weak  for  such  a  small  tumor.  In  these 
images,  we  use  e  =  0.1M  and  the  searching  range  [—14,  14]. 
Fig.  4(a)  corresponds  to  DAS-C,  where  the  tumor  is  buried  by 
interference  and  noise.  In  Fig.  4(b),  DAS -El  fails  to  detect  the 
tumor.  In  Fig.  4(c),  for  DAS-E2,  a  shadow  of  the  tumor  can  be 
seen.  In  Fig.  4(d),  for  RCB-E2,  most  of  the  clutters  are  cleared 
up  but  a  strong  clutter  shows  up  near  the  chest  wall.  Fig.  4(e)- 
4(h)  shows  the  results  of  peak  searching;  none  of  them  have 
false  tumors,  which  may  be  attributed  to  proper  corrections  of 
phase  aberrations.  Images  produced  by  ARMOR-P  in  Fig.  4(f) 
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Fig.  4.  Reconstructed  images  based  on  the  2-D  simulated  breast  model  with 
one  2-mm-diameter  tumor,  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB- 
E2,  with  s  —  0.1  AT.  (e)  DAS-P.  (f)  ARMOR-P,  with  e  =  0.1M.  (g)  DAS-PP. 
(h)  ARMOR-PP,  with  e  =  0.1M. 


Fig.  5.  Reconstructed  images  based  on  the  2-D  simulated  breast  model  with 
one  large  tumor  (1  cm  in  diameter).  The  white  circle  in  the  image  corresponds 
to  the  actual  shape  of  the  tumor,  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB- 
E2,  with  s  —  0.1  AT.  (e)  DAS-P.  (f)  ARMOR-P,  with  e  =  0.1M.  (g)  DAS-PP. 
(h)  ARMOR-PP,  with  e  =  0.1  M. 


and  by  ARMOR-PP  in  Fig.  4(h)  have  lower  sidelobe  levels  and 
higher  resolutions,  and  the  latter  has  a  higher  contrast  than  the 
former,  due  to  the  latter  using  the  peak-to-peak  difference  as  the 
intensity  measure. 

Fig.  5  shows  the  imaging  results  for  the  one  large  tumor 
(1  cm  diameter)  case.  Here,  we  set  e  =  0.1  M  and  the  search¬ 
ing  range  [—20,  20].  (Note  that  different  tumor  sizes  and  loca¬ 
tions  will  result  in  different  sound  speed  variations  in  the  breast 
model.)  The  white  circle  in  the  image  corresponds  to  the  ac¬ 
tual  contour  of  the  tumor.  Although  all  the  methods  can  detect 
the  tumor,  only  ARMOR  can  be  used  to  form  an  image  of  the 
tumor  with  the  best  agreement  with  the  actual  tumor  size  and 
location. 

By  plotting  a  map  (maps  are  not  shown  here  due  to  limited 
space)  of  the  values  of  p  used  in  ARMOR,  for  each  focal  point, 
we  find  that  at  the  tumor  locations,  fi  usually  takes  smaller  values 
than  that  at  other  locations. 


B.  Experimental  Results 

We  have  also  tested  ARMOR  and  DAS  on  two  sets  of  TAT 
experimental  data  from  mastectomy  specimens  [4]  obtained  by 
the  Optical  Imaging  Laboratory  at  the  Texas  A&M  University. 

The  two  data  sets  were  acquired  from  mastectomy  specimens 
using  a  TAT  system.  Microwave  sources  were  used  to  heat  the 
specimens  transiently.  In  the  experiment,  the  breast  specimen 
was  formed  to  a  cylindrical  shape  inside  a  plastic  bowl.  The 
bowl  was  immersed  in  ultrasound  coupling  medium  in  a  con¬ 
tainer.  For  breast  specimen  I,  the  acoustic  signals  were  recorded 
at  240  equally  spaced  scanning  stops  on  a  circular  track  of  radius 
12.9  cm.  The  thickness  of  this  specimen  was  about  4  cm  in  a 
round  plastic  bowl  of  17  cm  in  diameter.  This  lesion  was  diag¬ 
nosed  as  an  invasive  metaplastic  carcinoma  with  chondroid  and 
squamous  metaplasia.  The  size  of  the  tumor  was  measured  to  be 
35  mm  in  diameter  by  TAT,  and  36  mm  in  diameter  by  radiog¬ 
raphy  (see  [4]  for  details).  For  breast  specimen  II,  the  scanning 
radius  was  9.7  cm,  with  160  scanning  stops.  This  specimen  was 
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G) 


Fig.  6.  Reconstructed  images  for  breast  specimen  I.  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB-E2,  with  £  =  0.5 M.  (e)  DAS-R  (f)  ARMOR-P,  with 
e  =  0.5M.  (g)  DAS-PP.  (h)  ARMOR-PP,  with  e  =  0.5 M.  (i)  X-ray  image,  (j)  Inverse  solution. 


9  cm  thick  in  a  round  plastic  bowl  of  1 1  cm  in  diameter.  The 
lesion  in  the  specimen  was  diagnosed  as  infiltrating  lobular  car¬ 
cinoma;  the  size  of  the  tumor  was  about  20  mm  x  12  mm  on 
TAT  image,  and  about  26  mm  x  15  mm  on  the  radiography 
(see  [4]  for  more  details). 

First,  we  study  the  delay  difference  for  both  the  breast  speci¬ 
mens  to  get  a  qualitative  guide  for  choosing  the  searching  range 
in  Step  II  of  ARMOR.  The  results  are  shown  in  Fig.  3(c)  and 
3(d),  respectively.  Note  that  breast  specimen  II  has  a  larger  vari¬ 
ance  in  delay  differences  than  breast  specimen  I.  In  Fig.  3(c), 
70%  of  the  delay  differences  are  roughly  between  —23  to  23 
samples,  whereas  in  Fig.  3(d),  70%  of  the  delay  differences 
are  between  —40  and  40  samples.  Therefore  we  should  set  a 
larger  searching  range  for  breast  specimen  II  than  for  breast 
specimen  I. 


Fig.  6  shows  the  reconstructed  images  for  breast  specimen  I. 
In  the  following  images,  the  searching  range  was  set  to  [—3,  3] 
after  adjustment,  and  e  —  0.5M  for  all  the  RCBs  used  herein. 
In  Fig.  6(a),  for  DAS-C,  the  dark  region  shows  a  blurred  object 
corresponding  to  the  breast  tumor.  In  Fig.  6(b),  for  DAS-E1,  the 
light  region  shows  a  vague  boundary  of  the  tumor.  Fig.  6(c), 
for  DAS-E2,  and  6(d),  for  RCB-E2,  have  similar  performances. 
In  Fig.  6(e),  for  DAS-P,  and  6(f),  for  ARMOR-P,  a  dark  region 
with  a  clear  cut  has  a  good  correspondence  with  the  location 
and  shape  of  the  tumor  in  the  radiograph  [4].  In  Fig.  6(g),  for 
DAS-PP,  and  6(h),  for  ARMOR-PP,  not  only  a  clear  image  of  the 
tumor  is  obtained,  but  also  the  detailed  boundary  is  revealed.  For 
comparison,  the  images  from  X-ray  mammography,  considered 
the  “gold  standard”  of  breast  imaging,  and  the  exact  inverse 
solution  of  TAT  (see  [4]  for  more  details)  are  shown  in  Fig.  6(i) 
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Fig.  7.  Reconstructed  images  for  breast  specimen  II.  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB-E2,  with  £  =  0.5 M.  (e)  DAS-P.  (f)  ARMOR-P,  with 
£  =  0.5 M .  (g)  DAS-PP.  (h)  ARMOR-PP,  with  £  =  0.5 M .  (i)  X-ray  image,  (j)  Inverse  solution. 


and  6(j),  respectively.  We  give  Fig.  6  and  the  following  Fig.  7  in 
gray  scale  to  have  a  better  comparison  with  the  X-ray  images. 

Fig.  7  shows  the  reconstructed  images  for  breast  specimen 
II.  The  tumor  size  here  is  smaller,  and  a  high  level  of  interfer¬ 
ence  and  noise  is  present  in  the  recorded  data.  The  searching 
interval  is  eventually  adjusted  to  [—120,  120]  and  RCB  pa¬ 
rameter  6  =  0.5M.  In  Fig.  7(a),  for  DAS-C,  the  true  tumor  is 
barely  identifiable  from  the  surrounding  clutters.  The  DAS -El 
shown  in  Fig.  7(b)  and  the  DAS-E2  shown  in  Fig.  7(c)  provide 
higher  imaging  contrast  than  DAS-C  but  show  strong  clutter.  In 
Fig.  7(d),  for  RCB-E2,  a  false  tumor  shows  up,  which  demon¬ 
strates  the  need  for  robustness  in  the  presence  of  relatively  strong 
phase  distortion.  DAS-P  is  shown  in  Fig.  7(e)  and  ARMOR-P  is 
shown  in  Fig.  7(f).  DAS-PP  and  ARMOR-PP  produce  the  best 
images  in  Figs.  7(g)  and  7(h),  respectively,  with  the  location  and 


shape  of  the  tumor  consistent  with  the  radiograph  in  Fig.  7(i)  [4] . 
If  we  define  the  signal-to-background  ratio  (SBR)  (i.e.,  squaring 
the  pixel  values  of  the  image,  the  ratio  of  the  maximum  to  the 
total  sum  of  the  squared  values)  as  an  image  quality  measure¬ 
ment  metric,  ARMOR-PP  has  an  SBR  twice  that  of  DAS-PP, 
which  means  a  3  dB  gain  for  ARMOR-PP.  For  comparison,  the 
image  formed  by  the  exact  inverse  solution  of  TAT  (see  [4]  for 
more  details)  is  shown  in  Fig.  7(j). 

The  effects  of  the  uncertainty  parameter  6  in  ARMOR  is 
studied  in  our  next  example.  We  vary  5  of  RCB  used  in  ARMOR. 
The  imaging  results  for  breast  specimen  I,  shown  in  Fig.  8, 
are  consistent  with  our  previous  analysis:  when  5  is  large,  the 
performance  of  RCB,  in  Fig.  8(a),  is  close  to  that  of  DAS  in 
Fig.  6(g).  When  the  parameter  6  is  small,  as  shown  in  Fig.  8(c), 
the  resolution  is  improved  at  the  cost  of  robustness. 
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Fig.  8.  Effects  of  uncertainty  parameter  £  on  ARMOR-PP,  with  a  searching 
range  [-3,  3].  (a)  e  =  0.7 M.  (b)  £  =  0.5M.  (c)  £  =  0.3M. 
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Fig.  9.  Effects  of  the  searching  range  on  the  DAS-PP  images,  (a)  Searching 
range  [—20,  20].  (b)  Searching  range:  [—40, 40].  (c)  Searching  range  [—60,  60]. 
(d)  Searching  range:  [—80,  80]. 

In  our  last  example,  the  effect  of  the  searching-range  width  on 
the  imaging  quality  is  considered.  We  use  DAS-PP  as  an  exam¬ 
ple  since  it  shows  more  dependence  on  the  searching  range.  The 
conclusion  drawn  for  DAS  applies  to  ARMOR.  A  symmetric 
searching  range  centered  around  the  calculated  arrival  time  is 


used.  From  the  discussions  in  Section  IV,  we  know  that  there 
is  a  tradeoff  in  choosing  the  searching  range.  Clearly,  when  the 
searching  range  is  too  small,  such  as  in  Fig.  9(a),  we  miss  the 
true  peaks.  With  an  increase  in  the  searching  range,  the  image 
quality  becomes  gradually  better,  as  shown  in  Fig.  9(b)  and  9(c). 
However,  when  the  searching  range  passes  a  certain  threshold, 
with  too  much  interference  coming  into  the  searching  range,  the 
image  quality  degrades  because  of  increased  clutters,  as  shown 
in  Fig.  9(d). 

From  our  numerical  examples,  we  conclude  that  ARMOR  has 
higher  resolution  and  better  interference  rejection  capability  and 
more  robustness  against  wavefront  distortion  than  DAS.  Also, 
we  find  that  the  amplitude-based  measures  reveal  more  details  of 
the  tumor  in  the  reconstructed  images  than  their  energy-based 
counterparts.  The  energy-based  measures  are  not  sensitive  to 
phase  distortions;  however,  they  tend  to  blur  the  reconstructed 
images,  causing  loss  of  details  with  a  low-pass  filtering-like 
effect. 


VII.  Conclusion 

ARMOR  has  been  proposed  for  thermoacoustic  tomogra¬ 
phy.  ARMOR  is  robust  to  the  amplitude  and  phase  distortions 
in  the  recorded  signals  caused  by  the  acoustic  heterogeneity 
of  biological  tissues.  ARMOR  consists  of  three  steps:  in  the 
first  step,  ARMOR  uses  the  data-adaptive  robust  Capon  beam¬ 
forming  (RCB)  for  waveform  estimation;  in  the  second  step  of 
ARMOR,  a  simple,  yet  effective,  peak  searching  method  is  used 
to  mitigate  the  phase  distortion  in  the  estimated  waveform;  in 
the  third  step,  the  response  intensity  is  calculated  for  the  focal 
point  using  various  approaches,  among  which  the  peak-to-peak 
difference  measure  further  enhances  the  image  contrast.  Exam¬ 
ples  based  on  a  numerically  simulated  2-D  breast  model  and  two 
sets  of  experimentally  measured  data  from  human  mastectomy 
specimens  demonstrate  the  excellent  performance  of  ARMOR: 
high-resolution,  low  sidelobe  level,  and  much  improved  inter¬ 
ference  suppression  capability. 

Appendix 

Thermal  Acoustic  Simulations 

We  consider  the  microwave-induced  thermal  acoustic  simu¬ 
lation  in  two  steps.  In  the  first  step,  the  electromagnetic  field 
inside  the  breast  model  is  simulated  and  the  specific  absorp¬ 
tion  rate  (SAR)  distribution  is  calculated  based  on  the  simu¬ 
lated  electromagnetic  field.  The  second  step  is  for  the  acoustic 
wave  simulation,  where  the  SAR  distribution  obtained  in  the 
first  step  is  used  as  the  acoustic  pressure  source  through  the 
thermal  expansion  coefficient.  In  both  steps,  the  finite-difference 
time-domain  (FDTD)  method  [28]  is  used  for  the  simulation  ex¬ 
amples. 

The  2-D  electromagnetic  breast  model  used  is  as  shown  in 
Fig.  2(a).  A  narrow  electromagnetic  pulse  is  used  to  irradiate 
the  breast  from  the  top  of  the  model.  The  electromagnetic  field 
is  simulated  using  the  FDTD  method.  The  grid-cell  size  used  by 
FDTD  is  0.5  mm  x  0.5  mm  and  the  computational  region  is  ter¬ 
minated  by  perfectly  matched  layer  (PML)  absorbing  boundary 
conditions  [29]. 
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TABLE  II 

Nominal  Dielectric  Properties  of  Breast  Tissues  [24] 


Tissues 

Dielectric  Properties 

Permittivity  (F/m) 

Conductivity  (S/m) 

Immersion  Liquid 

9 

0 

Chest  Wall 

50 

7 

Skin 

36 

4 

Fatty  Breast  Tissue 

9 

0.4 

Nipple 

45 

5 

Glandular  Tissue 

11-15 

0.4-0. 5 

Tumor 

50 

4 

TABLE  III 

Acoustic  Parameters  for  Biological  Tissues 


Tissue 

p  (kg/m3) 

c  (m/s) 

a*  (dB/cm) 

(3  (l/°  C) 

Cp  (J/(°  C  •  kg)) 

Breast 

1020 

1510 

0.75/1-5 

3E-4 

3550 

Skin 

1100 

1537 

3.5 

3E-4 

3500 

Muscle 

1041 

1580 

0.57/ 

3E-4 

3510 

Tumor 

1041 

1580 

0.57/ 

3E-4 

3510 

*/  is  the  acoustic  frequency,  and  the  unit  is  in  megahertz. 


The  SAR  distribution  is  given  as  [30] 


SAR(r) 


a(r)E2  (r) 

2p(r) 


(26) 


where  cr(r)  is  the  conductivity  of  the  biological  tissues  at  loca¬ 
tion  r,  E( r)  is  the  electric  field  at  location  r,  and  p( r)  is  the 
mass  density  of  the  biological  tissues  at  location  r. 

In  the  microwave-induced  TAI  system,  the  microwave  energy 
is  small,  and  as  a  result,  the  acoustic  pressure  field  induced  by 
the  microwave  is  also  small.  So,  the  nonlinear  acoustic  effect 
does  not  need  to  be  considered  in  the  TAI  system.  The  two  basic 
linear  acoustic  wave  generation  equations  are  [9] 

p2u(r,  t)  =  —Vp(r,  t )  (27) 

and 

\  d  d 

V  •  u(r ,t)  = - ^—p(r,t)+ap(r,t)+f3—T(r,t)  (28) 

pcz  ot  ot 

where  u(r,  t)  is  the  acoustic  velocity  vector,  p(r,  t)  is  the  acous¬ 
tic  pressure  field,  p  is  the  mass  density,  a  is  the  attenuation 
coefficient,  (3  is  the  thermal  expansion  coefficient,  and  T(r,  t) 
is  the  temperature.  The  values  for  these  acoustic  properties  for 
different  breast  tissues  are  listed  in  Table  III  [25]. 

Because  the  duration  of  the  microwave  pulse  is  much  shorter 
than  the  thermal  diffusion  time,  thermal  diffusion  can  be  ne¬ 
glected  [9],  and  the  thermal  equation  is 

CpEr(r,t)  =  SAR(r,t)  (29) 

where  Cp  is  the  specific  heat.  Substituting  (29)  into  (28)  gives 


\  d  8 

V  •  u(r,  t)  = - 2  w.p(r,  t)  +  ap( r,  t)  +  -LsAR(r,  t). 

pc  Ol  Up 

(30) 

FDTD  is  used  again  to  compute  the  thermal  acoustic  wave  based 
on  (27)  and  (30). 

The  breast  model  for  the  acoustic  simulation  is  shown  in 
Fig.  2(b),  which  is  constructed  similarly  to  the  model  for  elec¬ 
tromagnetic  simulation.  An  acoustic  sensor  array  deployed  uni¬ 


formly  around  the  breast  model  is  used  to  record  the  thermal 
acoustic  signals.  The  grid-cell  size  used  by  the  acoustic  FDTD 
is  0.1  mm  x  0.1  mm  and  the  computational  region  is  termi¬ 
nated  by  PML-absorbing  boundary  conditions.  Note  that  the 
size  of  the  FDTD  cell  for  the  acoustic  simulation  is  much  finer 
than  that  of  the  FDTD  cell  for  the  electromagnetic  simulation 
because  the  wavelength  of  an  acoustic  wave  is  much  smaller 
than  that  of  a  microwave.  The  SAR  distribution  data  is  interpo¬ 
lated  to  achieve  a  desired  grid  resolution  for  the  acoustic  breast 
model. 
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Waveform  Diversity  Based  Ultrasound  System  for 
Hyperthermia  Treatment  of  Breast  Cancer 

Bin  Guo,  Member,  IEEE ,  and  Jian  Li,  Fellow,  IEEE 


A  bstract-  In  this  letter,  we  present  a  new  waveform-diver¬ 
sity-based  ultrasound  hyperthermia  technique  for  the  treatment 
of  breast  cancer.  Waveform  diversity  offers  a  new  paradigm  for 
beampattern  design.  By  choosing  a  proper  covariance  matrix  of 
the  transmitted  waveforms  under  the  uniform  elemental  power 
constraint,  the  ultrasound  system  can  provide  a  focal  spot  matched 
to  the  entire  tumor  region,  and  meanwhile,  minimize  the  impact 
to  the  surrounding  healthy  breast  tissues.  As  shown  in  our  2-D 
numerical  simulations,  this  method  has  better  acoustic  power  de¬ 
position  than  the  existing  methods,  and  can  provide  the  necessary 
temperature  gradients  required  for  the  effective  hyperthermia 
treatment  of  breast  cancer. 

Index  Terms—  Beampattern  design,  breast  cancer,  finite-differ¬ 
ence  time  domain  (FDTD),  ultrasound  hyperthermia,  waveform 
diversity. 


I.  Introduction 

BREAST  CANCER  is  the  most  common  nonskin  malig¬ 
nancy  in  women  and  the  second  leading  cause  of  female 
cancer  mortality  [1].  There  are  over  200000  new  cases  of  in¬ 
vasive  breast  cancer  diagnosed  each  year  in  the  United  States, 
and  one  out  of  every  seven  women  in  the  United  States  will  be 
diagnosed  with  breast  cancer  in  their  life  time.1 

The  development  of  breast  cancer  imaging  techniques,  such 
as  microwave  imaging  [2],  [3],  ultrasound  imaging  [4],  [5], 
thermal  acoustic  imaging  [6],  and  magnetic  resonance  imaging 
(MRI),  has  improved  the  ability  to  visualize  and  accurately 
locate  the  breast  tumor  without  the  need  for  surgery  [7].  This 
has  led  to  the  probability  of  noninvasive  local  hyperthermia 
treatment  of  breast  cancer.  Many  studies  have  been  performed 
to  demonstrate  the  effectiveness  of  the  local  hyperthermia  on 
the  treatment  of  breast  cancer  [8],  [9].  A  challenge  in  the  local 
hyperthermia  treatment  of  breast  cancer  is  to  heat  the  malignant 
tumors  to  a  temperature  above  43  °C  for  about  30-60  min,  but 
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to  maintain  a  low-temperature  level  in  the  surrounding  healthy 
breast  tissue  region. 

There  are  two  major  classes  of  local  hyperthermia  tech¬ 
niques:  microwave  hyperthermia  [10]  and  ultrasound  hy¬ 
perthermia  [11].  The  penetration  of  microwave  is  poor  in 
biological  tissues.  Moreover,  the  focal  spot  generated  by  mi¬ 
crowave  is  undesirable  at  the  normal/cancerous  tissues  interface 
because  of  the  long  wavelength  of  the  microwave.  Ultrasound 
can  achieve  much  better  penetration  depths  than  microwave. 
However,  because  the  acoustic  wavelength  is  very  short,  the 
focal  spot  generated  by  ultrasound  is  very  small  (millimeter  or 
submillimeter  in  diameter)  compared  to  the  large  tumor  region 
(centimeter  in  diameter  on  average).  Thus,  many  focal  spots 
are  required  for  complete  tumor  coverage,  and  this  results  in  a 
long  treatment  time  and  missed  cancer  cells. 

In  this  letter,  we  present  a  waveform-diversity-based  ul¬ 
trasound  hyperthermia  technique  for  the  treatment  of  breast 
cancer.  Waveform  diversity  is  a  new  beampattern  design  tech¬ 
nique  recently  proposed  for  multiple-input-multiple-output 
(MIMO)  radar  (see  [12]  and  the  references  therein).  Unlike 
the  standard  phased-array  technique,  transmitting  multiple 
different  waveforms  via  its  transducers  offers  more  flexibility 
for  transmit  beampattern  design.  By  designing  the  transmitted 
signal  cross-correlation  matrix  under  the  uniform  elemental 
power  constraint,  the  waveform  diversity  can  be  exploited  to 
maximize  the  power  deposition  at  the  entire  tumor  region  while 
minimizing  the  impact  on  the  surrounding  healthy  tissue  region. 

To  validate  our  algorithm,  we  develop  a  2-D  breast  model 
with  an  embedded  tumor.  The  model  includes  the  breast  tissue, 
skin,  and  chest  wall.  The  finite-difference  time-domain  (FDTD) 
method  is  used  to  simulate  the  acoustic  field  and  the  temperature 
distribution  within  the  breast.  We  show  with  numerical  simula¬ 
tions  that  the  proposed  method  can  provide  the  necessary  tem¬ 
perature  gradients  required  for  the  effective  hyperthermia  treat¬ 
ment  of  the  tumor  and  maintain  a  low-temperature  level  at  the 
surrounding  healthy  tissue  region. 


II.  Waveform-Diversity-Based  Ultrasound 
Hyperthermia 

We  consider  an  ultrasound  hyperthermia  system  as  shown  in 
Fig.  1.  Let  1*0  denote  the  center  location  of  the  tumor,  which  is 
assumed  to  be  estimated  accurately  a  priori  using  breast  cancer 
imaging  techniques.  There  are  M  acoustic  transducers  deployed 
around  the  breast  at  locations  rm  (ra  =  1,  2,  . . . ,  M).  Let 
Xm(n)  (n  =  1,  2,  . . . ,  TV)  denote  the  discrete-time  baseband 
signal  transmitted  by  the  rath  acoustic  transducer,  where  TV  de¬ 
notes  the  number  of  samples  of  each  transmitted  signal  pulse. 
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We  assume  that  the  transmitted  acoustic  signals  are  narrow- 
band  and  each  acoustic  transducer  is  omnidirectional.  The  base¬ 
band  signal  at  a  location  r  inside  the  breast  can  be  described  as 


e-j27T/0Tm(r) 

y(r,n)  =  \  J - rn -7oXm(n),  n=l,2,...,N 

^Ikm-rll172 

(1) 

where  /o  is  the  carrier  frequency 

Tm(  r)  =  -  (2) 

c 


is  the  time  needed  by  the  signal  emitted  via  the  rath  transducer 
to  arrive  at  the  location  r,  with  c  being  the  sound  speed  inside 
the  breast  tissues,  and  l/||rm  —  rH1/2  is  the  propagation  atten¬ 
uation  of  the  acoustic  wave.  Let 


a(r) 


"  ei27r/0ri(r)  ej27r/oT2(r) 

||ri  —  rH1/2  ||r2  —  r\\1/2 


ei27r/0rM(r) 

1 1 I'M  -  I'll1/2 


(3) 


be  the  steering  vector  with  (-)T  denoting  the  transpose,  and  let 
x(n)  =  [xi{n)  x2(n)  ■■■  xM(n)]T .  (4) 

Equation  (1)  can  be  rewritten  as 

y(r,n)  =  a*(r)x(n),  n  =  l,2,...,N  (5) 

where  (•)*  denotes  the  conjugate  transpose. 

The  power  of  the  transmitted  signals  at  location  r,  which  is 
also  called  the  transmit  beampattern  [12],  is  given  by 

P(  r)  =  E{y(r,n)y*(r,n)}  =  a*(r)Ra(r)  (6) 

where  R  is  the  covariance  matrix  of  x(n),  i.e., 

R  =  E{x(n)x*(n)}.  (7) 

The  transmit  beampattern  is  a  function  of  the  location  r. 

The  purpose  of  our  waveform-diversity  technique  is  to  focus 
the  acoustic  power  onto  the  entire  tumor  region  while  mini¬ 
mizing  the  peak  power  level  at  the  surrounding  healthy  breast 


tissue  region.  The  corresponding  beampattern  design  problem  is 
to  choose  the  covariance  matrix  R  under  the  uniform  elemental 
power  constraint 

Rmm  =  jj,  ra  =  1,  2,  . . . ,  M  (8) 

where  denotes  the  (ra,  ra)th  element  of  R,  and  b  is  the 
total  transmitted  power,  to  achieve  the  following  goals: 

1)  achieve  a  predetermind  main-beam  width  matching  the 
entire  tumor  region  (be  within  10%  of  the  power  de¬ 
posited  at  the  tumor  center); 

2)  minimize  the  peak  sidelobe  level  in  a  prescribed  region 
(the  surrounding  healthy  breast  tissue  region). 

This  problem  can  be  formulated  as 


min  —t 

i,R 

s.t.  a*(r0)Ra(r0)  —  a*(/^)Ra(/i)  >  t  V/t  G  Hs 
a*(i/)Ra(i/)  >  0.9a*(ro)Ra(ro)  V  v  e  Hr 

a*(i/)Ra(i/)  <  l.la*(r0)Ra(r0)  V  v  €  fir 


R  >  0 

Rmm=Yj,  ra  =  1,  2,  . . . ,  M  (9) 


where  Hr  and  H b  denote  the  tumor  region  and  the  surrounding 
healthy  breast  tissue  regions  (sidelobe  region),  respectively. 

As  shown  in  [12],  this  beampattern  design  problem  is  a 
semidefinite  program  (SDP)  and  can  be  efficiently  solved  in 
polynomial  time  using  public  domain  software.  Once  R  is  de¬ 
termined,  a  signal  sequence  (x(n) }  that  has  R  as  its  covariance 
matrix  can  be  synthesized  as 


x(n)  =  R1//2^(n),  n  =  1,  2,  . . . ,  N  (10) 


where  {w(n)}  is  a  sequence  of  independent  identically  dis¬ 
tributed  (i.i.d.)  random  vectors  with  mean  zero  and  covariance 
matrix  I,  and  R1/2  denotes  a  square  root  of  R. 

By  transmitting  x(n)  given  in  (10)  using  the  acoustic  trans¬ 
ducer  array,  we  can  approximately  get  a  desired  high  acoustic 
power  deposition  matching  the  entire  tumor  region  while  min¬ 
imizing  the  power  deposition  at  the  surrounding  healthy  breast 
tissue  region. 


III.  Model  and  Numerical  Results 

A.  Breast  Model  and  Simulation 

For  simulation  purposes,  a  2-D  breast  model  is  established, 
as  shown  in  Fig.  1.  The  breast  model  is  a  10-cm-diameter  semi¬ 
circle,  which  includes  breast  tissues,  skin,  and  chest  wall.  The 
acoustic  properties  of  the  breast  tissues  within  the  breast  are 
assumed  random  with  a  variation  of  ±5%  around  the  nominal 
values.  A  16-mm-diameter  tumor  is  embedded  below  the  skin 
with  the  tumor  center  location  being  (x  =  0  mm  and  y  = 
50  mm).  There  are  51  acoustic  transducers  deployed  uniformly 
around  the  breast  model.  The  distance  between  the  neighboring 
acoustic  transducers  is  1.5  mm  (half  wavelength  of  the  carrier 
frequency).  Acoustic  wave  with  frequency  500  kHz  is  used  as 
the  carrier  frequency. 
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Fig.  2.  Beampattern.  (a)  Waveform  diversity,  (b)  Phased  array. 


TABLE  I  TABLE  II 

Typical  Acoustic  Properties  of  Breast  Tissues  Typical  Thermal  Properties  of  Breast  Tissues 


Tissues 

p  (kg/m6) 

c  (m/s) 

a  (dB/cm) 

Breast  Tissue 

1020 

1510 

0.26 

Skin 

1100 

1537 

1.0 

Chest  Wall 

1041 

1580 

0.28 

Tumor 

1041 

1580 

0.28 

Tissues 

K  (  w  ) 

a  (W_\ 

c  (  J  \ 

K  t m-°C ' 

A  lm3  J 

B  UC.m3  J 

C  t kq-°C ) 

Breast  Tissue 

0.499 

480 

2700 

3550 

Skin 

0.376 

1620 

9100 

3500 

Chest  Wall 

0.564 

480 

2700 

3510 

Tumor 

0.564 

480 

2700 

3510 

The  two  basic  linear  acoustic  wave  generation  equations  are  tissues.  The  thermal  model  is  based  on  the  bioheat  equation 
[13],  [14]  nm 


pj^ufr,  ^  =  _  0  (H) 

and 

1  d 

V  •  u(r, t)  = - ^—p(r,t)  +  ap{r,t)  (12) 

pcz  at 

where  u(r ,t)  is  the  acoustic  velocity  vector,  p(r,t)  is  the 
acoustic  pressure  field,  p  is  the  mass  density,  and  a  is  the 
attenuation  coefficient.  The  nominal  values  for  these  acoustic 
properties  for  different  breast  tissues  are  listed  in  Table  I  [4], 
[15],  [16].  The  values  for  the  tumor  are  approximated  using 
those  for  muscle  because  we  cannot  find  the  values  specific 
to  the  tumor.  FDTD  is  used  to  compute  the  acoustic  field 
distribution  based  on  (11)  and  (12).  More  details  about  FDTD 
for  acoustic  simulations  can  be  found  in  [17]  and  [18]. 

Once  the  acoustic  pressure  is  calculated,  the  acoustic  power 
deposition  at  location  r,  denoted  as  Q( r),  is  given  as  [14] 


Q(r)  =  -b(r)|2.  (13) 

pc 

After  obtaining  the  acoustic  power  deposition,  the  2-D 
thermal  model,  corresponding  to  the  2-D  acoustic  models, 
is  used  to  calculate  the  temperature  distribution  in  the  breast 


\(K(r)X7T}±A(r)+Q(r)-B(r)(T-TB)  =  C(r)p(r)  —  (14) 

where  K(r)  is  the  thermal  conductivity,  A(r)  is  metabolic  heat 
production,  B(r)  represents  the  heat  exchange  mechanism  due 
to  capillary  blood  perfusion,  C{ r)  is  the  specific  heat,  and  Tb 
is  the  blood  temperature,  which  can  be  assumed  as  the  body 
temperature.  The  thermal  properties  for  our  breast  model  are 
listed  in  Table  II.  More  detailed  discussions  can  be  found  in  [19] . 

The  thermal  models  are  also  simulated  using  the  FDTD 
method  [20].  The  body  temperature  and  the  environmental 
temperature  are  set  at  36.8  °C  and  20  °C,  respectively.  The 
convective  boundary  condition  is  used  at  the  skin  surface. 


B.  Numerical  Results 


We  demonstrate  the  performance  of  our  waveform-diver¬ 
sity-based  method  via  several  numerical  examples.  For  compar¬ 
ison  purposes,  the  conventional  delay-and-sum  (DAS)-based 
phased-array  beamforming  method  is  also  applied  to  the  same 
model,  and  its  results  are  compared  with  those  of  our  wave- 
form-diversity-based  method. 

DAS-based  phased-array  beamformer  transmits  the  same 
waveform  using  a  weight  vector 


w  = 


afro) 

|a(r0)||2' 


(15) 
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Fig.  3.  Power  deposition,  (a)  Waveform  diversity,  (b)  Phased  array. 


C 


Fig.  4.  Temperature  distribution,  (a)  Waveform  diversity,  (b)  Phased  array. 


The  corresponding  beampattern  is 

P(r)  =  |a*(r)  w|2.  (16) 

Fig.  2  shows  the  calculated  beampattern  within  the  breast 
model.  The  blue  lines  in  these  figures  mark  the  boundary  of  the 
breast  and  the  tumor.  Fig.  2(a)  is  the  beampattern  due  to  our 
waveform  diversity  technique  which  is  calculated  by  using  (6) 
with  the  optimal  covariance  matrix  R  determined  by  using  (9). 
The  figure  shows  that  the  3-dB  main  beam  is  matched  to  the 
tumor  region  well,  and  the  sidelobe  level  is  low.  Fig.  2(b)  is  the 
DAS  beampattern  which  is  calculated  using  (16).  It  is  shown 
that  the  DAS  beampattern  is  very  narrow,  and  only  focuses  at 
the  center  region  of  the  tumor. 

Fig.  3(a)  and  (b)  shows  the  acoustic  power  densities  within 
the  breast  model  for  the  waveform  diversity  technique  and  DAS, 
respectively.  The  dots  mark  the  locations  of  the  acoustic  trans¬ 
ducers.  It  is  shown  that  the  acoustic  power  densities  in  Fig.  3 


agree  with  the  beampatterns  in  Fig.  2  very  well,  and  our  wave¬ 
form  diversity  technique  gives  a  focal  spot  matching  the  entire 
tumor  region. 

Fig.  4  shows  the  temperature  distributions  within  the  breast 
model.  Fig.  4(a)  is  the  result  of  waveform  diversity,  which  shows 
that  the  entire  tumor  region  is  heated  to  a  temperature  greater 
than  43  °C  while  maintaining  the  surrounding  healthy  tissues  at 
a  low-temperature  level  (below  40  °C).  As  a  comparison,  DAS 
only  heats  a  very  small  region  to  a  temperature  greater  than 
43  °C  at  the  center  of  the  tumor. 

IV.  Conclusion 

In  this  letter,  we  have  presented  a  new  waveform-diversity- 
based  ultrasound  hyperthermia  technique  for  the  treatment  of 
breast  cancer.  By  choosing  the  covariance  matrix  of  the  trans¬ 
mitted  waveforms  properly,  this  method  can  provide  a  focal  spot 
matching  the  entire  tumor  region  while  minimizing  the  impact 
on  the  surrounding  healthy  breast  region.  As  shown  with  2-D 
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simulation  examples,  this  method  has  better  acoustic  power  de¬ 
position  than  its  conventional  counterpart,  and  can  provide  the 
necessary  temperature  gradients  required  for  the  effective  hy¬ 
perthermia  of  breast  cancer. 
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Abstract-  Microwave-induced  thermal  acoustic  imaging  (TAI) 
is  a  promising  early  breast  cancer  detection  technique,  which  com¬ 
bines  the  advantages  of  microwave  stimulation  and  ultrasound 
imaging  and  offers  a  high  imaging  contrast,  as  well  as  high  spatial 
resolution  at  the  same  time.  A  new  multifrequency  microwave-in¬ 
duced  thermal  acoustic  imaging  scheme  for  early  breast  cancer 
detection  is  proposed  in  this  paper.  Significantly  more  information 
about  the  human  breast  can  be  gathered  using  multiple  frequency 
microwave  stimulation.  A  multifrequency  adaptive  and  robust 
technique  (MART)  is  presented  for  image  formation.  Due  to  its 
data-adaptive  nature,  MART  can  achieve  better  resolution  and 
better  interference  rejection  capability  than  its  data-independent 
counterparts,  such  as  the  delay-and-sum  method.  The  effective¬ 
ness  of  this  procedure  is  shown  by  several  numerical  examples 
based  on  2-D  breast  models.  The  finite-difference  time-domain 
method  is  used  to  simulate  the  electromagnetic  field  distribution, 
the  absorbed  microwave  energy  density,  and  the  thermal  acoustic 
field  in  the  breast  model. 

Index  Terms—  Breast  cancer  detection,  finite-difference  time-do- 
main  (FDTD)  methods,  multifrequency  adaptive  and  robust 
technology  (MART),  robust  capon  beamforming  (RCB),  thermal 
acoustic  imaging  (TAI). 


I.  Introduction 

BREAST  cancer  is  the  most  common  nonskin  malignancy 
in  women  and  the  second  leading  cause  of  female  cancer 
mortality  [1].  There  are  over  200000  new  cases  of  invasive 
breast  cancer  diagnosed  each  year  in  the  U.  S.,  and  one  out  of 
every  seven  women  in  the  U.S.  will  be  diagnosed  with  breast 
cancer  in  their  life  time  (the  American  Cancer  Society,  2006) 
and  early  diagnosis  is  key  to  surviving  breast  cancer  [2].  Mi¬ 
crowave  imaging  is  a  method  for  early  breast  cancer  detection 
[3]-[9],  which  exploits  the  significant  contrast  in  dielectric 
properties  between  normal  and  cancerous  tissues  [10]— [12]. 
However,  it  is  difficult  for  microwave  imaging  techniques  to 
achieve  good  (submillimeter)  spatial  resolution  because  of 
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Fig.  1.  Model  of  microwave-induced  TAI  for  breast  cancer  detection. 


its  long  wavelength  [13].  Ultrasound  is  another  option  which 
offers  a  high  spatial  resolution  because  of  its  short  acoustic 
wavelength  [14]— [16] .  However,  the  contrast  in  acoustic  prop¬ 
erties  between  normal  and  tumor  tissues  is  very  small  due  to 
both  being  soft  tissues. 

Microwave-induced  thermal  acoustic  imaging  (TAI)  com¬ 
bines  the  advantages  of  microwave  stimulation  and  ultrasound 
imaging  [13],  which  offers  a  high  imaging  contrast  (due  to  the 
significantly  different  dielectric  properties  of  tumor  and  normal 
breast  tissues),  as  well  as  high  spatial  resolution  (due  to  the  low 
propagation  velocity  or  the  short  wavelength  of  acoustic  waves 
in  biological  tissues)  at  the  same  time.  To  use  microwave-in¬ 
duced  TAI  techniques  for  breast  cancer  imaging,  a  microwave 
source  with  a  short  duration  time  is  used  to  irradiate  the  breast, 
as  shown  in  Fig.  1.  The  normal  breast  tissues,  as  well  as  tumors, 
absorb  microwave  energy  and  emanate  thermal  acoustic  waves 
by  thermoelastic  expansion.  It  is  well-known  that  malignant 
breast  tissue  has  a  higher  water  content  [1],  [10],  [12],  [17], 
with  a  much  higher  conductivity  than  normal  breast  tissues 
(with  low  water  content).  As  a  result,  the  microwave  energy  ab¬ 
sorbed  by  tumor  and  normal  breast  tissues  will  be  significantly 
different  and  a  stronger  acoustic  wave  will  be  produced  by  the 
tumor.  The  acoustic  waves  generated  in  this  manner  carry  the 
information  about  the  microwave  energy  absorption  properties 
of  the  tissues  under  irradiation.  The  thermal  acoustic  waves 
propagate  out  of  the  breast  and  are  recorded  by  an  acoustic 
sensor  array  placed  around  the  breast.  The  tumor  locations  can 
be  accurately  determined  since  the  received  acoustic  signals 
from  the  malignant  tumors  are  at  higher  levels,  hence  aiding 
image  construction. 

During  the  last  decade,  several  research  groups  have  been 
working  on  the  microwave-induced  TAI  of  biological  tissues 
[1 8]— [23] .  The  microwave  frequency  used  ranges  from  400  MHz 
[22]  to  3  GHz  [13].  Image  reconstruction  algorithms  include  the 
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widely  used  delay-and-sum  (DAS)  method  [20],  [23],  the  fre¬ 
quency-domain  inverse  method  [24],  [25],  and  the  time-domain 
inverse  method  [13],  [19]. 

Microwave-induced  TAI  does,  however,  present  several 
challenges.  First,  the  human  breast  is  large  in  size,  usually 
has  an  irregular  shape  if  not  compressed,  and  is  covered  with 
a  2-mm-thick  skin  with  dielectric  properties  significantly 
different  from  the  normal  breast  tissues.  Moreover  the  breast 
tissue  is  far  from  homogeneous  because  the  dielectric  proper¬ 
ties  of  glandular  tissue  are  different  from  that  of  breast  fatty 
tissue.  All  these  factors  make  it  difficult  to  approximate  the 
back  propagation  properties  of  thermal  acoustic  signals  inside 
the  breast.  Due  to  the  slow  acoustic  wave  propagation  speed  or 
short  wavelength  in  biological  tissues,  the  errors  on  the  order 
of  millimeters  in  determining  the  acoustic  signal  propagation 
path  lengths  will  severely  degrade  the  image  quality. 

In  this  paper,  a  multifrequency  microwave-induced  TAI 
system  is  proposed  which  remedy  the  problems  mentioned 
above.  Instead  of  using  a  single  frequency  microwave  source, 
as  generally  done  by  other  research  groups  in  this  field,  here 
a  multiple  frequency  source  is  used,  since  the  desired  thermal 
acoustic  signals  can  be  induced  by  microwave  sources  operating 
at  a  wide  range  of  frequencies.  We  show  in  this  paper  that  the 
rich  information  collected  from  the  multifrequency  stimulation 
can  help  mitigate  the  challenges  mentioned.  The  multifre¬ 
quency  microwave-induced  thermal  acoustic  signals  will  offer 
higher  signal-to-noise  ratio  (SNR)  and  higher  imaging  contrast 
than  single-frequency  microwave-induced  thermal  acoustic  sig¬ 
nals  since  much  more  information  about  the  human  breast  can 
be  harvested  from  the  multiple  stimulating  frequencies  within 
the  microwave  frequency  band.  Furthermore,  the  interference 
due  to  inhomogeneous  breast  tissue  can  be  suppressed  more 
effectively  when  multifrequency  microwave-induced  thermal 
acoustic  signals  are  used  for  image  reconstruction  since  more 
information  about  breast  tissue  can  be  used  by  the  adaptive 
image  reconstruction  algorithms. 

Another  challenge  encountered  by  microwave-induced  TAI 
is  the  need  to  develop  accurate  and  robust  image  reconstruction 
methods.  DAS  is  a  widely  used  reconstruction  algorithm  in 
medical  imaging.  This  method  is  data-independent  and  tends  to 
suffer  from  poor  resolution  and  high  sidelobe  level  problems. 
Data-adaptive  approaches,  such  as  the  recently  introduced 
robust  Capon  beamforming  (RCB)  [26],  [27]  method,  can  have 
much  better  resolution  and  much  better  interference  rejection 
capabilities  than  their  data-independent  counterpart.  Several 
medical  imaging  algorithms  [4],  [5],  [28],  [29]  based  on  RCB 
have  been  developed  and  used  for  microwave  imaging  and 
thermal  acoustic  imaging.  Good  performances  of  these  algo¬ 
rithms  have  been  reported. 

We  present  a  multifrequency  adaptive  and  robust  technique 
(MART)  based  on  RCB  for  multifrequency  microwave-induced 
TAI.  There  are  three  stages  in  our  MART.  In  Stage  I,  RCB  is 
used  to  estimate  the  thermal  acoustic  responses  from  the  grid 
locations  within  the  breast  for  each  stimulating  microwave 
frequency.  Then,  in  Stage  II,  a  scalar  acoustic  waveform  at 
each  grid  location  is  estimated  based  on  the  response  estimates 
for  all  stimulating  frequencies  from  Stage  I.  Finally,  in  Stage 
III,  the  positive  peak  and  the  negative  peak  of  the  estimated 


acoustic  waveform  at  each  grid  location  are  determined,  and 
the  peak-to-peak  difference  is  computed  and  referred  to  as  the 
image  intensity. 

To  validate  the  effectiveness  of  the  proposed  algorithm,  we 
develop  a  2-D  inhomogeneous  breast  model,  which  includes 
skin,  breast  fatty  tissues,  glandular  tissues,  and  the  chest  wall. 
Small  tumors  are  set  below  the  skin.  The  finite-difference  time- 
domain  (FDTD)  method  is  used  to  simulate  the  electromagnetic 
field  inside  the  breast  tissues  [30],  [31].  The  specific  absorp¬ 
tion  rate  (SAR)  distribution  is  calculated  based  on  the  simu¬ 
lated  electromagnetic  field  [32],  [33].  Then  FDTD  is  used  again 
to  simulate  the  propagation  of  the  microwave-induced  thermal 
acoustic  waves  [34],  [35]. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  II,  the  microwave  frequency  properties  of  human 
breast  are  described.  A  proper  microwave  frequency  band  for 
multifrequency  microwave-induced  TAI  is  also  given  in  this 
section.  MART  is  proposed  for  image  formation  in  Section  III. 
In  Section  IV,  2-D  electromagnetic  and  acoustic  breast  models 
are  developed.  The  electromagnetic  and  acoustic  simulation 
methods  are  also  presented  in  this  section.  Imaging  results 
based  on  numerical  examples  are  provided  in  Section  V. 
Section  VI  concludes  the  paper. 

II.  Microwave  Properties  of  Human  Breast 

A.  Cutoff  Frequency  of  Human  Breast 

In  a  microwave-induced  TAI  system,  the  biological  tissues 
should  be  heated  by  microwave  sources  in  a  uniform  manner, 
otherwise  thermal  acoustic  signals  will  be  induced  by  a  nonuni¬ 
form  microwave  energy  distribution,  resulting  in  images  diffi¬ 
cult  to  interrupt.  It  is  well-known  that  high-order  electromag¬ 
netic  field  modes  will  be  excited  in  a  media  if  the  microwave 
works  at  a  frequency  higher  than  a  cutoff  frequency  of  the  media 
[36],  and  the  microwave  energy  distribution  is  nonuniform  at 
high-order  modes  [37] .  To  minimize  the  nonuniform  microwave 
energy  distribution  inside  the  breast  caused  by  the  high-order 
electromagnetic  modes,  the  microwave  source  should  work  at  a 
frequency  below  a  certain  cutoff  frequency. 

To  estimate  the  cutoff  frequency  for  the  human  breast,  we 
consider  the  simplified  breast  model  shown  in  Fig.  2(a)  con¬ 
sisting  of  a  semicircular  dielectric  waveguide  with  a  perfect 
magnetic  conductor  (PMC)  at  the  bottom  of  the  semicircle.  Re¬ 
call  that  the  tangential  components  of  the  magnetic  field  are  zero 
on  the  surface  at  the  PMC.  The  PMC  assumption  is  reasonable 
because  the  permittivity  of  the  chest  wall  (er  =  50)  is  much 
greater  than  that  of  the  normal  breast  tissues  (er  =  9).  In  cir¬ 
cular  dielectric  waveguide,  if  an  electromagnetic  mode  has  a 
field  distribution  whose  tangential  magnetic  field  components 
are  zero  at  the  center  line  of  the  circular  waveguide,  as  shown 
in  Fig.  2(b),  the  introduction  of  a  PMC  at  the  center  line  of  the 
circular  waveguide  will  not  significantly  change  the  boundary 
conditions  and,  hence,  will  not  significantly  alter  the  mode  dis¬ 
tribution.  The  modes  in  the  semicircular  dielectric  waveguide 
can,  thus,  be  estimated  by  determining  the  modes  in  a  corre¬ 
sponding  circular  dielectric  waveguide. 
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Fig.  2.  Simplified  breast  model,  (a)  Semicircular  dielectric  waveguide  with 
PEC,  and  (b)  corresponding  circular  dielectric  waveguide. 


The  dominant  mode  of  a  circular  dielectric  waveguide  is 
the  HE  11  mode,  the  cutoff  frequency  of  which  is  zero.  The 
electromagnetic  field  distribution  is  near  uniform  at  this  mode. 
The  dominant  mode  is  followed  by  the  TE01,  TM01,  and 
HE21  modes.  These  modes  are  degenerate,  and  have  a  cutoff 
frequency  given  by  [36] 


fc  = 


__XoiCo__ 
27 xa\Jer  —  1 


(1) 


where  Co  is  the  speed  of  light  in  free  space,  xoi  —  2.405  is  the 
first  root  of  the  Bessel  function  of  the  first  kind  of  order  zero 
(Jo(xoi)  =  0),  a  and  er  are  the  radius  and  average  permit¬ 
tivity  of  the  circular  dielectric  waveguide,  respectively.  TM01 
and  HE21,  as  well  as  the  interference  between  them,  satisfy  the 
zero  tangential  magnetic  field  component  condition  at  the  center 
line  of  the  circular  waveguide.  These  modes  can  also  exist  in  the 
semicircular  dielectric  waveguide.  By  substituting  the  parame¬ 
ters  of  the  breast  model  into  (1),  we  obtain  the  cutoff  frequency 
of  the  semicircular  breast  model  to  be 


fc  = 


2.405  •  3  x  108 
2tt  •  0.05  • 


=  812  MHz 


(2) 


where  we  have  used  a  —  5  cm  and  er  =  9  as  typical  values 
for  human  breast.  Consequently,  the  stimulating  microwave  fre¬ 
quency  for  the  TAI  system  should  be  below  812  MHz. 


B.  Microwave  Energy  Absorption  Properties  of  Breast  Tissues 
and  Tumor 

It  is  well-known  that  the  complex  relative  dielectric  properties 
of  a  medium  can  be  expressed  as 


^ r  ^ r 


TABLE  I 

Cole-Cole  Parameters  for  Biological  Tissues 


Tissue 

Breast 

Skin 

Muscle 

Tumor 

£oo 

2.5 

4.0 

4.0 

4.0 

a 

0.01 

0.0002 

0.2 

0.7 

Aeri 

3.0 

32.0 

50.0 

50.0 

n  (ps) 

17.68 

7.23 

7.23 

7.0 

OL\ 

0.1 

0.0 

0.1 

0.0 

As2 

15 

1100 

7000 

0 

T2  (ns) 

63.66 

32.48 

353.68 

N/A 

OL2 

0.1 

0.2 

0.1 

N/A 

As3 

5.0E4 

0 

1.2E6 

0 

T3  (MS) 

454.7 

N/A 

318.31 

N/A 

G3 

0.1 

n/a 

0.1 

N/A 

A^4 

2.0E7 

n/a 

2.5E7 

0 

T4  (ms) 

13.26 

n/a 

2.274 

N/A 

OL4. 

0.0 

n/a 

0.0 

N/A 

where  e'r  is  the  relative  permittivity  and  e"  is  the  out-of-phase 
loss  factor  which  can  be  written  as 


with  a  being  the  total  conductivity,  £q  the  free  space  permit¬ 
tivity,  and  u )  the  electromagnetic  frequency.  The  tissue  absorp¬ 
tion  property  of  the  electromagnetic  wave  energy  is 

Q{v)=1-<j\E{v)\2  (5) 

which  is  a  function  of  the  total  conductivity  and  the  electric  field 
inside  the  tissue.  If  we  assume  that  the  microwave  energy  distri¬ 
bution  is  uniform  inside  the  breast  in  a  TAI  system,  the  absorp¬ 
tion  of  the  microwave  energy  by  the  breast  is  characterized  by 
the  total  conductivity  of  the  breast  tissues 


=  s^Equj. 


(6) 


Hence,  instead  of  using  the  attenuation  coefficient  a ,  as  used  in 
[23],  in  this  paper,  we  study  the  absorption  properties  of  breast 
tissues  using  the  total  conductivity  a. 

The  dielectric  properties  of  biological  tissues  can  be  accu¬ 
rately  modeled  by  the  Cole-Cole  equation  [38] 


K 

£r(id)  =  Sqo  +  ^ 

i= 1 


A  £j 

l  +  (juny-^ 


+ 


^0 

jue  o 


(7) 


where  K  is  the  order  of  the  Cole-Cole  model,  -&m  is  the 
high-frequency  permittivity,  t7;  is  the  relaxation  time,  A Ei  is 
the  pole  amplitude,  ai  (0  <  a.i  <  1)  is  a  measure  of  the  broad¬ 
ening  of  dispersion,  and  ao  is  the  static  ionic  conductivity.  The 
Cole-Cole  parameters  for  skin,  breast  fatty  tissue,  chest  wall 
(mainly  consisting  of  muscle),  as  well  as  tumor  are  listed  in 
Table  I  [39],  [40].  Because  we  cannot  find  the  values  specific  to 
the  tumor,  the  dielectric  properties  of  the  tumor  is  approximated 
using  a  Debye  model  [3],  [41],  which  is  a  special  case  of  the 
Cole-Cole  model. 

Substituting  (7)  into  (6),  we  obtain  the  total  conductivity  of 
the  breast  tissue  as  follows: 


t(lo)  =  -imag  + 


K 


A  £j 


(3) 


jue  o 


(8) 
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Fig.  3.  Total  conductivity  of  normal  breast  tissues  and  tumor  as  a  function  of 
frequency. 


Fig.  4.  Ratio  of  conductivity  between  tumor  and  normal  breast  tissue  as  a  func¬ 
tion  of  frequency. 


which  is  a  function  of  the  stimulating  microwave  frequency, 
where  imag(-)  denotes  the  imaginary  part  of  the  complex 
relative  permittivity.  Fig.  3  gives  the  total  conductivity  of 
breast  fatty  tissue  and  tumor  over  a  frequency  band  from 
100-1000  MHz.  Note  that  the  total  conductivity  increases  with 
the  microwave  frequency,  which  means  that  more  microwave 
energy  is  absorbed  and  converted  to  heat  by  tissues  at  higher 
microwave  frequency  region,  or  in  other  words,  the  SNR  is 
higher  in  the  received  thermal  acoustic  signals  at  higher  stim¬ 
ulating  microwave  frequency  region.  On  the  other  hand,  the 
penetration  at  higher  microwave  frequencies  is  smaller  because 
the  tissues  are  lossy.  We  define  the  conductivity  ratio  between 
the  tumor  and  the  normal  breast  tissue  as 


r<r(u)  = 


^tumor(^) 
®  breast  (^) 


(9) 


and  plot  it  as  a  function  of  frequency  in  Fig.  4.  A  high  con¬ 
ductivity  ratio  means  that  more  microwave  energy  is  absorbed 
and  converted  to  heat  by  tumor  than  by  normal  breast  tissues. 
In  other  words,  the  higher  the  conductivity  ratio,  the  higher 
the  imaging  contrast.  Fig.  4  shows  that  the  imaging  contrast 
is  higher  at  the  lower  microwave  frequency  region  because  the 
conductivity  ratio  decreases  with  the  frequency. 

These  microwave  energy  absorption  properties  of  breast 
tissues  and  tumor  motivate  us  to  consider  inducing  thermal 
acoustic  signals  with  different  microwave  frequencies.  By 


taking  into  account  the  aforementioned  cutoff  frequency  given 
in  (2),  we  choose  a  frequency  range  from  200-800  MHz.  The 
frequency  step  is  100  MHz,  with  a  total  of  seven  frequencies. 
Note  that  wideband  antenna  techniques  should  be  used  for 
the  practical  implementation  because  the  frequency  range  is 
wide.  However,  since  the  exciting  microwave  frequency  is 
stepped,  an  antenna  with  a  broad  instantaneous  bandwidth 
is  not  required.  Another  advantage  of  using  multiple  fre¬ 
quencies  for  stimulation  is  that  more  information  about  the 
inhomogeneous  breast  tissues  will  be  harvested  from  the  mul¬ 
tifrequency  microwave-induced  thermal  acoustic  signals.  The 
microwave  energy  distribution  inside  the  breast  model  is  not 
uniform  because  the  human  breast  is  inhomogeneous  media, 
and  thermal  acoustic  signals  will  be  induced  by  the  inhomo¬ 
geneous  energy  distribution.  These  thermal  acoustic  signals 
will  appear  as  clutter  in  the  resulting  images.  However,  the 
inhomogeneous  microwave  energy  distributions  are  different  at 
different  stimulating  frequencies  because  of  the  different  mi¬ 
crowave  wavelengths  in  breast  tissues.  When  a  multifrequency 
microwave  source  is  used  for  TAI,  the  thermal  acoustic  clutter 
induced  by  the  inhomogeneous  breast  tissues  can  be  suppressed 
by  our  adaptive  and  robust  imaging  algorithm. 

III.  Multifrequency  Adaptive  and  Robust  Technology 
(MART)  for  Breast  Cancer  Imaging 

We  consider  a  multifrequency  microwave-induced  TAI 
system  as  shown  in  Fig.  1,  where  an  acoustic  sensor  array  is 
arranged  on  a  semicircle  relatively  close  to  the  breast  skin.  The 
location  of  each  acoustic  sensor  is  r j  (j  =  1,  •  •  • ,  TV),  where 
N  is  the  number  of  acoustic  sensors.  Assume  that  M  —  7  mi¬ 
crowave  sources  with  different  frequencies  are  used  to  irradiate 
the  breast  model.  Let  Pij(t)  (i  =  1,  •  •  • ,  M;  j  =  1,  •  •  • ,  N\ 
t  =  0,  •  •  • ,  T  —  1)  denote  the  thermal  acoustic  signal  induced 
by  the  ith  frequency  and  received  by  the  jth  acoustic  sensor, 
where  T  is  the  recording  time  which  is  sufficiently  long  to 
allow  acoustic  sensors  to  record  all  responses  from  the  breast. 
Our  goal  is  to  detect  the  tumor  by  reconstructing  an  image  of 
the  thermal  acoustic  response  intensity  J(r)  as  a  function  of 
scan  location  r  within  the  breast. 

A.  Data  Preprocessing 

Because  the  breast  skin,  breast  tissues,  chest  wall,  and  tumor 
absorb  the  microwave  energy  and  convert  the  energy  to  heat,  all 
of  them  produce  thermal  acoustic  signals.  The  received  thermal 
acoustic  waveforms  include  the  responses  from  the  tumor,  as 
well  as  from  other  healthy  breast  tissues.  The  thermal  acoustic 
signals  generated  by  the  skin  are  much  stronger  than  those  by  a 
small  tumor  because  of  the  high  conductivity  of  the  skin  and  the 
acoustic  sensors  being  very  close  to  the  skin.  We  must  remove 
the  skin  responses  to  enhance  the  tumor  responses.  Because  the 
distances  between  the  acoustic  sensors  and  the  nearest  breast 
skin  are  similar  to  one  another,  the  signals  recorded  by  var¬ 
ious  sensors  have  similar  skin  responses.  Hence,  we  can  remove 
the  skin  response  by  subtracting  out  a  fixed  calibration  signal 
from  all  received  signals.  This  calibration  signal  can  be  obtained 
simply  by  averaging  the  recorded  signals  from  all  channels. 

Let  Xij(t)  denote  the  signals  after  subtracting  out  the  calibra¬ 
tion  signal.  To  process  the  signals  coherently  for  a  focal  point 
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at  r,  we  align  the  signals  Xij(t)  by  time  shifting  each  signal  a 
number  of  samples  rij( r).  The  discrete  time  delay  between  r 
and  the  jth  acoustic  sensor  can  be  calculated  as 


nA r)  = 


(10) 


where  |_7j  stands  for  rounding  to  the  greatest  integer  less  than 
7,  ||i*j  —  r||  is  the  distance  between  17  and  r,  c  is  the  velocity 
of  the  acoustic  wave  propagating  in  breast  tissues,  and  At  is  the 
sampling  interval,  which  is  assumed  to  be  sufficiently  small.  The 
time- shifted  signals  are  denoted  as 


r)  =  Xij  ( t  +  nj( r)) ,  t  =  -rij( r),  ■  ■  ■  ,T  -  nj( r). 

(11) 

After  time  shifting,  the  acoustic  signals  from  the  imaging  lo¬ 
cation  r  are  aligned  so  that  they  all  start  approximately  from 
time  t  =  0  for  all  channels.  Now  the  aligned  signals  are  win¬ 
dowed  by 


Window(Z)  = 


1,  0  <1  <  L  —  1 
0,  otherwise 


(12) 


to  isolate  the  signals  from  the  focal  point  at  r.  The  windowed 
signals  are  denoted  as  Xij(L,  r),  /  =  0,  •••,  £  —  1,  where  LAt  is 
the  approximate  duration  of  the  thermal  acoustic  pulse,  which 
can  be  determined  from  the  pulse  duration  of  the  pulsed  mi¬ 
crowave  source. 

Attenuation  exists  when  acoustic  waves  propagate  within  the 
breast.  This  attenuation  has  two  parts:  the  attenuation  due  to  the 
lossy  media  and  the  propagation  attenuation.  Thus,  the  attenu¬ 
ation  of  the  tumor  responses  at  various  channels  are  different 
because  of  the  different  distances  between  the  imaging  position 
r  and  the  acoustic  sensors.  For  the  2-D  case  considered  here,  the 
compensation  factor  for  the  jth  channel  is  given  by 


Fig.  5.  Data  cube  model.  In  Stage  I,  MART  slices  the  data  cube  for  each  fre¬ 
quency  index.  RCB  is  applied  to  each  data  slice  to  estimate  the  corresponding 
waveform. 


MART  is  a  three-stage  time-domain  signal  processing  algo¬ 
rithm.  In  Stage  I,  MART  slices  the  data  cube  corresponding  to 
each  frequency  index,  and  processes  each  data  slice  by  the  RCB 
to  obtain  the  thermal  acoustic  waveform  estimate  for  each  stim¬ 
ulating  frequency.  Then,  in  Stage  II,  a  scalar  waveform  is  es¬ 
timated  from  all  frequencies  based  on  the  waveform  estimates 
from  Stage  I.  Finally,  the  positive  peak  and  the  negative  peak 
of  the  estimated  thermal  acoustic  waveform  from  Stage  II  are 
found  in  Stage  III.  The  peak-to-peak  difference  is  calculated  as 
the  image  intensity  at  the  focal  point  at  r.  The  details  of  all  three 
stages  are  given  below. 

1 )  Stage  I:  In  Stage  I,  MART  approximates  the  data  model 
as 


yi(Z)  =  a iSi{l)  +  e;(Z)  (16) 


Kj(r)  =  exp  (a\\rj  -  r||)  •  ||rj  -  r||1/2  (13) 

where  the  first  term  of  the  right  side  of  (13)  compensates  for  the 
attenuation  due  to  the  lossy  media,  and  the  second  term  com¬ 
pensates  for  the  geometric  attenuation.  The  compensated  signal 
can  be  calculated  as 

yij(l,r)  =  Kj(r)  ■  Xij(l,r),  1  =  0,  (14) 


B.  Multifrequency  Adaptive  and  Robust  Technology  (MART) 

Without  loss  of  generality,  we  consider  imaging  at  a  generic 
location  r  only.  For  notational  convenience,  we  drop  the  depen¬ 
dence  of  yij  (/,  r)  on  r,  and  simply  denote  it  as  Now  we 

consider  the  data  model 

Vi.j  (0  —  $i,j  (0  &i.  j  (0  (1^) 

where  represents  the  tumor  response  and  e^y(Z)  rep¬ 

resents  the  residual  term,  which  includes  the  noise  and  inter¬ 
ference  from  breast  skin,  chest  wall,  and  other  responses.  The 
structure  of  the  data  model  is  a  data  cube  as  shown  in  Fig.  5. 


where  y;(Z)  =  ■  ■  ■ ,  2/i,jv(Z)]T  and  e;(Z)  = 

[<%;l(Z),  ‘ ' '  >  e*,w(Z)]T-  The  scalar  waveform  st(l)  denotes 
the  thermal  acoustic  signal  generated  at  the  focal  location  r 
corresponding  to  the  7th  stimulating  frequency.  The  vector  is 
referred  to  as  the  array  steering  vector,  which  is  approximately 
equal  to  1atxi  =  [1,  •  •  • ,  1]T  since  all  the  signals  have  been 
aligned  temporally  and  their  attenuation  compensated  for  in  the 
preprocessing  step.  The  residual  e^(Z)  is  the  noise  and  interfer¬ 
ence  term,  which  is  assumed  uncorrelated  with  the  signal. 

There  are  two  assumptions  made  to  write  the  model  given  in 
(16).  First,  the  steering  vector  is  assumed  to  vary  with  the  mi¬ 
crowave  frequency  ( i )  but  nearly  constant  with  the  time  sample 
(/).  Second,  we  assume  that  the  thermal  acoustic  signal  wave¬ 
form  depends  only  on  the  microwave  frequency  ( i )  but  not  on 
the  acoustic  sensor  (  j ) .  The  truth,  however,  is  that  the  steering 
vector  is  not  exactly  known  as  it  changes  slightly  with  both  the 
stimulating  frequency  and  time  due  to  array  calibration  errors 
and  other  factors.  The  signal  waveform  can  also  vary  slightly 
with  both  the  stimulating  frequency  and  acoustic  sensor,  due  to 
the  inhomogeneous  and  frequency-dependent  medium  within 
the  breast.  The  two  aforementioned  assumptions  simplify  the 
problem  slightly.  They  cause  little  performance  degradations 
when  used  with  our  adaptive  and  robust  algorithm. 
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In  practice,  the  true  steering  vector  in  (16)  is  not  ljvxi-  We 
assume  that  the  true  steering  vector  a i  lies  in  the  vicinity  of  the 
assumed  steering  vector  a  =  1  jvxi ,  and  the  only  knowledge  we 
have  about  a i  is  that 


|a;-a||2  <*,  (17) 

where  e\  is  a  user  parameter,  which  may  be  determined  de¬ 
pending  on  the  various  errors  discussed  previously. 

The  true  steering  vector  a^  can  be  estimated  via  the  following 
covariance  fitting  approach  of  RCB  [26],  [27] 


ma xSf  subject  to  Ry,  —  >  0 

%  5  1 

||a*  -  a||2  <  6!  (18) 

where  Sf  is  the  power  of  the  signal  .s.,(Z)  and 

R-y  ^  !  y^yd/ly/V)  (19) 

J  1=0 

is  the  sample  covariance  matrix.  The  above  optimization 
problem  can  be  solved  as  described  in  [26],  and  the  estimated 
true  steering  vector  is  denoted  here  as  a i . 

To  obtain  the  signal  waveform  estimate,  we  pass  the  received 
signals  through  a  Capon  beamformer  [27],  [42].  The  weight 
vector  of  the  beamformer  is  determined  by  using  the  estimated 
steering  vector  a^  in  the  following  expression: 


R T1 


w  a  = 


Y 


a^R 


-C 


(20) 


Then  the  estimated  signal  waveform  corresponding  to  the  ith 
stimulating  frequency  is 


k(l)  =  w/y;(7).  (21) 


By  repeating  the  aforementioned  process  for  i  =  1  through 
i  =  M,  we  obtain  the  complete  set  of  M  waveform  estimates 

s(Z)  =  [s1(/),---,sm(0]T.  (22) 


2)  Stage  II:  Since  the  stimulating  microwave  sources  with 
various  frequencies  are  assumed  to  have  the  same  power,  we  as¬ 
sume  that  the  thermal  acoustic  responses  from  the  tumor  at  dif¬ 
ferent  stimulating  frequencies  have  nearly  identical  waveforms. 
Note  that  the  thermal  acoustic  responses  induced  by  the  inho¬ 
mogeneous  microwave  energy  distribution  (due  to  the  inhomo¬ 
geneous  breast  tissues)  are  different  at  different  stimulating  fre¬ 
quencies.  This  means  that  the  elements  of  the  vector  s (Z)  are  all 
approximately  equal  to  an  unknown  scalar  signal  s(Z),  and  the 
noise  and  interference  term  can  be  assumed  uncorrelated  with 
this  signal.  In  Stage  II  of  MART,  we  adopt  the  data  model 


where  as  is  approximately  equal  to  Imxi-  However,  the 
“steering  vector”  may  again  be  imprecise,  and,  hence,  RCB  is 
needed  again. 

As  we  did  in  Stage  I,  we  assume  that  the  only  knowledge 
about  as  is  that 


as||2  <  e2, 


(24) 


where  as  =  1m  xi  is  the  assumed  steering  vector,  and  62  is  a 
user  parameter.  Again,  the  true  steering  vector  as  can  be  esti¬ 
mated  via  the  covariance  fitting  approach 

max#2  subject  to  Rs  —  £2asaJ  >  0 

62  ,as 

||as  —  a,||2  <  e2  (25) 


where  82  is  the  power  of  the  signal  s(Z),  and 

r  »  =  2  j;s(osT(o 

■  1=0 


(26) 


is  the  sample  covariance  matrix. 

After  obtaining  the  estimated  steering  vector  as,  we  obtain 
the  adaptive  weight  vector  and  the  estimated  signal  waveform, 
respectively,  as 


ajRUas 


(27) 


and 


s(t )  =  wTs  (t).  (28) 

3)  Stage  III:  Because  the  thermal  acoustic  pulse  is  usually 
bipolar:  a  positive  peak,  corresponding  to  the  compression 
pulse,  and  a  negative  peak,  corresponding  to  the  rarefaction 
pulse  [43],  we  use  the  peak-to-peak  difference  as  the  response 
intensity  for  the  imaging  location  r  in  the  third  stage  of  MART. 
Compared  with  other  energy  or  amplitude  based  response 
intensity  estimation  methods,  peak-to-peak  difference  can  be 
used  to  improve  imaging  quality  with  little  additional  compu¬ 
tation  costs. 

The  positive  and  negative  peak  values  of  the  estimated  wave¬ 
form  for  the  focal  location  r  will  be  searched  based  on  the 
estimated  waveform  (28)  obtained  in  Stage  II.  Because  of  the 
nonuniform  sound  speed  in  biological  tissues,  the  arrival  time 
of  the  acoustic  pulse  generated  at  location  r  cannot  be  calcu¬ 
lated  accurately.  However,  it  was  reported  in  [18]  that,  when 
the  heterogeneity  is  weak,  such  as  in  breast  tissues,  amplitude 
distortion  caused  by  multipath  is  not  severe.  We  assume  that  the 
original  peak  remains  a  peak  in  the  estimated  waveform,  and  the 
positive  and  negative  peak  values  of  the  thermal  acoustic  pulse 
can  be  searched  as 

P+  =  max  <  max  |s(Z)}  ,  0  1  (29) 

( /G[A  i .  A-2,]  J 


s(Z)  =  a  ss(l)  +  es(l) 


(23) 


and 


P  =  min  <  min  {£(/)},  0 

(^G[Ai,A2] 


(30) 
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Fig.  6.  Breast  model  for  thermal  acoustic  simulation,  (a)  Model  for  electro¬ 
magnetic  simulation  and  (b)  model  for  acoustic  simulation. 


where  [Ai,  A2]  E  [0,1/]  is  the  searching  range.  Here  Ai  and 
A2  are  user  parameters,  and  the  details  on  how  to  choose  them 
can  be  found  in  [29]. 

After  the  positive  and  negative  peak  values  are  found,  the 
response  intensity  for  the  focal  point  at  location  r  is  given  as 

I(r)  =  P+  -P~.  (31) 


IV.  Modeling  and  Simulations 

We  consider  2-D  breast  models  simulated  in  two  steps.  In 
the  first  step,  the  electromagnetic  field  inside  the  breast  model 
is  simulated  and  the  SAR  distribution  is  calculated  based  on 
the  simulated  electromagnetic  field.  The  second  step  is  for  the 
acoustic  wave  simulation,  where  the  SAR  distribution  obtained 
in  the  first  step  is  used  as  the  acoustic  pressure  source  through 
the  thermal  expansion  coefficient.  In  both  steps,  the  FDTD 
method  is  used  for  the  simulation  examples. 

A.  Electromagnetic  Model  and  Simulation 

For  simulation  purposes,  the  2-D  electromagnetic  breast 
model  used  is  as  shown  in  Fig.  6(a).  The  breast  model  is  a 
10  cm  in  diameter  semicircle,  which  includes  the  skin,  breast 
fatty  tissue,  glandular  tissues,  and  chest  wall  (muscle).  A 
1-mm-diameter  tumor  is  embedded  below  the  skin.  The  di¬ 
electric  properties  of  the  breast  tissues  as  well  as  tumor  at  the 
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Fig.  7.  Gaussian  modulated  microwave  source. 


microwave  frequency  fi(i  =  1,  •  •  • ,  M)  were  calculated  based 
on  the  Cole-Cole  model  in  (7).  The  dielectric  properties  of  the 
normal  breast  fatty  tissue  are  assumed  random  with  a  variation 
of  ±10%  around  the  nominal  values.  The  dielectric  constants 
of  glandular  tissues  are  between  er  =  11  and  er  =  15. 

Fig.  7  shows  a  Gaussian  modulated  electromagnetic  wave 
used  to  irradiate  the  breast  from  the  top  of  the  model,  as  shown  in 
Fig.  6(a).  The  time  duration  for  the  Gaussian  pulse  is  1  ps.  The 
electromagnetic  field  is  simulated  using  the  FDTD  method  [30], 
[31].  The  grid  cell  size  used  by  FDTD  is  0.5  x  0.5  mm  and  the 
computational  region  is  terminated  by  perfectly  matched  layer 
(PML)  absorbing  boundary  conditions  [44],  [45]. 

The  SAR  distribution  is  given  as  [32],  [33] 


SAR(r) 


^(r)£2(r) 

2p{v) 


(32) 


where  cr(r)  is  the  conductivity  of  the  biological  tissues  at  loca¬ 
tion  r,  E( r)  is  the  electric  field  at  location  r,  and  p(r)  is  the 
mass  density  of  the  biological  tissues  at  location  r. 


B.  Acoustic  Model  and  Simulation 

In  the  microwave-induced  TAI  system,  the  microwave  energy 
is  small,  and  as  a  result,  the  acoustic  pressure  field  induced  by 
microwave  is  also  small.  So,  the  nonlinear  acoustic  effect  does 
not  need  to  be  considered  in  the  TAI  system.  For  example,  it  is 
shown  in  [46]  that  the  temperature  rise  is  about  0.1  mK  and  the 
acoustic  pressure  change  is  only  about  100  Pa  in  the  microwave- 
induced  TAI  system.  The  shock  distance  in  breast  tissues  is  [47] 

B  9  B  2 

v  —  2 A  t  P0°  \  _  2 A  m  P0C  t  C 

2(s2  +  1)  P  2(^:  +  1)  P  /max 

(33) 

where  (B/A)(&  10)  is  the  nonlinear  factor  of  the  breast  tis¬ 
sues,  po(~  1000  kg/m3)  is  the  mass  density  of  the  breast  tis¬ 
sues,  and  c(«  1500  m/s)  is  the  sound  speed  inside  the  breast 
tissues  [14].  p  is  the  acoustic  pressure  rise,  and  Amin  and  /max 
are  the  minimal  acoustic  wavelength  and  the  maximal  acoustic 
frequency  of  the  thermal  acoustic  signal,  respectively.  For  our 
breast  model,  the  acoustic  pressure  rise  is  p  —  100  Pa,  and  the 
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TABLE  II 

Acoustic  Parameters  for  Biological  Tissues. 

(*  /  Is  the  Acoustic  Frequency  and  the  Unit  Is  Megahertz) 


Tissue 

Breast 

Skin 

Muscle 

Tumor 

p  {kg/m6) 

1020 

1100 

1041 

1041 

c  {m/s) 

1510 

1537 

1580 

1580 

a *  {dB/cm) 

0.75/lb 

3.5 

0.57/ 

0.57/ 

P  (1/°C) 

3E-4 

3E-4 

3E-4 

3E-4 

Cp  ( J/(°C  ■  kg)) 

3550 

3500 

3510 

3510 

maximal  acoustic  frequency  is  /max  =  500  KHz.  By  substi¬ 
tuting  the  parameters  into  (33),  we  obtain  the  shock  distance  in 
breast  tissues  to  be 

_  2A  m  P0°2  m  C 
2(^  +  l)  V  /max 

5  1000  •  15002  1500 

”  2(5  +  1)  100  500  x  103 

=  2.8  x  104  m.  (34) 

Because  the  size  of  our  breast  model  is  only  10  cm,  which  is 
much  smaller  than  the  shock  distance,  it  is  reasonable  to  ig¬ 
nore  the  nonlinear  acoustic  effect  in  the  microwave-induced  TAI 
system. 

The  two  basic  linear  acoustic  wave  generation  equations  are 
[13] 

P^u(r,  t)  =  —  Vp(r,  t)  (35) 

and 

\  d  d 

V  •  u(r, i)  =  --^—p(r,t)  +  ap(r,t)+/3—T(r,t)  (36) 

where  u(r,t)  is  the  acoustic  velocity  vector,  p(r,t)  is  the 
acoustic  pressure  field,  p  is  the  mass  density,  a  is  the  at¬ 
tenuation  coefficient,  6  is  the  thermal  expansion  coefficient, 
and  T(r,  t)  is  the  temperature.  The  values  for  these  acoustic 
properties  for  different  breast  tissues  are  listed  in  Table  II  [14], 
[46],  [48]-[50].  The  attenuation  coefficient  is  calculated  with 
/  =  0.15  MHz.  The  values  for  the  tumor  are  approximated 
using  the  values  for  muscle  because  we  cannot  find  the  values 
specific  to  the  tumor. 

Because  the  duration  of  the  microwave  pulse  is  much  shorter 
than  the  thermal  diffusion  time,  thermal  diffusion  can  be  ne¬ 
glected  [13],  and  the  thermal  equation  is 

CpJ^T(r,  i)  =  SAR(r,  t)  (37) 

where  Cp  is  the  specific  heat.  Substituting  (37)  into  (36)  gives 

\  d  3 

V-u(r,£)  = - 2— p(r,t)  +  ap(r,t)  +  -^-SAR(r,t ).  (38) 

pc  ot  Up 

FDTD  is  used  again  to  compute  the  thermal  acoustic  wave 
based  on  (35)  and  (38).  More  details  about  FDTD  for  acoustic 
simulations  can  be  found  in  [34],  [35],  and  [5 1]— [57]. 

The  breast  model  for  the  acoustic  simulation  is  constructed 
similarly  to  the  model  for  electromagnetic  simulation.  The 


Thermal  Acoustic  Signals  from  Tumor 


Normalized  Spectrums  of  The  Thermal  Acoustic  Signals 


Fig.  8.  Thermal  acoustic  signals  at  different  stimulating  frequencies  /  =  200, 
400,  600,  and  800  MHz.  (a)  Thermal  acoustic  responses  from  tumor  only  and 
(b)  the  normalized  spectrums  of  the  signals  in  (a). 


velocities  of  the  normal  fatty  breast  tissue  are  also  assumed 
random  with  a  variation  of  ±5%  around  average  values,  as 
shown  in  Fig.  6(b).  An  acoustic  sensor  array  with  35  elements 
deployed  uniformly  around  the  breast  model  is  used  to  record 
the  thermal  acoustic  signals.  The  distance  between  neighboring 
acoustic  sensors  is  4  mm.  The  grid  cell  size  used  by  the  acoustic 
FDTD  is  0.1  x  0.1  mm  and  the  computational  region  is  ter¬ 
minated  by  PML  absorbing  boundary  conditions  [55]-[57]. 
Note  that  the  size  of  the  FDTD  cell  for  acoustic  simulation 
is  much  finer  than  that  of  the  FDTD  cell  for  electromagnetic 
simulation  because  the  wavelength  of  an  acoustic  wave  is  much 
smaller  than  that  of  a  microwave.  The  SAR  distribution  data 
is  interpolated  to  achieve  the  designed  grid  resolution  for  the 
acoustic  breast  model. 

V.  Numerical  Examples 

At  the  beginning  of  this  section,  the  typical  microwave-in¬ 
duced  thermal  acoustic  responses  from  the  tumor  are  plotted 
in  Fig.  8(a)  for  stimulating  frequencies  of  /  =  200,  400,  600, 
and  800  MHz.  The  signals  are  simulated  based  on  the  afore¬ 
mentioned  2-D  breast  model.  To  obtain  the  signals,  we  perform 
the  simulation  twice  at  each  stimulating  frequency,  with  and 
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Fig.  9.  Imaging  results  for  the  case  of  a  single  1  mm-diameter  tumor, 
(a)  MART,  (b)  DAS,  (c)  SART  at  stimulating  frequency  /  =  200  MHz, 
(d)  SART  at  stimulating  frequency  /  =  400  MHz,  (e)  SART  at  stimulating  fre¬ 
quency  /  =  600  MHz,  and  (f)  SART  at  stimulating  frequency  /  =  800  MHz. 


Fig.  10.  Imaging  results  for  the  two  1.5-mm-diameter  tumors  case,  (a)  Breast 
model,  (b)  MART,  (c)  DAS,  (d)  zoom  in  of  (b),  (e)  SART  at  stimulating  fre¬ 
quency  /  =  300  MHz,  and  (f)  SART  at  stimulating  frequency  /  =  700  MHz. 


without  the  tumor,  and  record  the  acoustic  signals  in  an  acoustic 
sensor.  The  difference  of  the  two  received  signals  is  referred  to 
as  the  thermal  acoustic  response  only  from  the  tumor  at  the  stim¬ 
ulating  frequency.  It  can  be  seen  that  the  thermal  acoustic  re¬ 
sponses  from  the  tumor  at  different  stimulating  frequencies  are 
similar  to  one  another.  The  figure  also  shows  that  the  thermal 
acoustic  signals  are  wideband  bipolar  pulses,  with  a  large  pos¬ 
itive  peak  and  a  large  negative  peak.  Fig.  8(b)  shows  the  nor¬ 
malized  spectra  of  the  acoustic  signals  corresponding  to  the  ex¬ 
citations  in  Fig.  8(a).  It  is  seen  that  the  frequency  range  of  the 
acoustic  signals  is  about  from  1  to  400  KHz.  The  dominant  band 
(3-dB  band)  of  the  signals  ranges  from  10  to  180  KHz,  and  the 
corresponding  acoustic  wavelength  ranges  from  150  to  8  mm  in 
the  breast  tissues. 

Several  numerical  examples  are  used  in  this  section  to  demon¬ 
strate  the  effectiveness  of  MART.  The  thermal  acoustic  sig¬ 
nals  are  simulated  based  on  the  2-D  breast  model  with  tumor 
only.  For  comparison  purposes,  the  DAS  method  is  applied  to 
the  same  data  set  also.  We  also  present  the  imaging  results  for 
the  single-frequency  microwave-induced  TAI  at  different  stim¬ 
ulating  frequencies.  The  corresponding  image  reconstruction 
method  is  referred  to  as  the  single-frequency  adaptive  and  robust 
technique  (SART),  which  is  similar  with  MART  but  without 
Stage  II  of  MART.  In  SART,  the  RCB  is  used  to  estimate  the 
thermal  acoustic  waveform  at  a  certain  stimulating  frequency 
just  like  in  Stage  I  of  MART.  Then  the  peak  search  method  used 
in  MART  Stage  III  is  applied  to  the  estimated  waveform  to  de¬ 
termine  the  image  intensity. 


In  the  first  example,  a  1 -mm-diameter  tumor  is  embedded  in 
the  breast  model  at  the  location  (x  =  70  mm,  y  =  60  mm). 
This  is  the  challenging  case  of  early  breast  cancer  detection 
because  of  the  small  tumor  size.  Fig.  9(a)  and  (b)  shows  the 
imaging  results  for  MART  and  DAS,  respectively.  The  tumor 
is  shown  clearly  in  the  MART  image  [Fig.  9(a)],  and  the  size 
and  location  of  the  tumor  is  accurate.  Because  of  the  high  side- 
lobe,  poor  resolution,  and  poor  interference  rejection  capability 
of  the  DAS  method,  the  tumor  is  essentially  missed  by  DAS  as 
shown  in  Fig.  9(b).  Fig.  9(c)-(f)  shows  the  imaging  results  for 
SART  at  the  stimulating  frequencies  /  =  200,  400,  600,  and 
800  MHz,  respectively.  The  figures  show  that  SART  can  de¬ 
termine  the  tumor  correctly,  but  some  clutter  shows  up  in  the 
SART  images.  Note  that  the  clutter  shows  up  at  different  loca¬ 
tions  with  different  stimulating  frequencies.  By  comparing  the 
images  for  MART  and  SART,  it  can  be  seen  that  the  clutter  are 
effectively  suppressed  by  MART  when  multiple  stimulating  fre¬ 
quencies  are  used. 

In  the  second  numerical  example,  two  small  1.5-mm-di- 
ameter  tumors  are  set  inside  the  breast  model  as  shown  in 
Fig.  10(a).  Their  locations  are  at  (x  =  70  mm,  y  =  60  mm)  and 
(x  =  75  mm,  y  =  62.5  mm).  The  distance  between  the  two 
tumors  is  4  mm.  The  imaging  results  using  MART  and  DAS 
are  shown  in  Fig.  10(b)  and  (c),  respectively.  The  two  tumors 
are  seen  clearly  in  the  MART  image.  To  show  them  clearly 
we  zoom  in  onto  the  tumor  locations  in  Fig.  10(d),  where  the 
two  black  circles  mark  the  actual  sizes  and  locations  of  the 
two  tumors.  It  is  shown  that  MART  can  be  used  to  determine 
the  locations  and  sizes  of  the  two  tumors  accurately.  The  DAS 
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image  contains  much  clutter.  The  two  tumors  cannot  be  sepa¬ 
rated  clearly  in  the  DAS  image  because  of  the  poor  resolution 
of  DAS.  Fig.  10(e)  and  (f)  shows  the  imaging  results  of  SART 
at  stimulation  frequencies  /  =  300  and  700  MHz,  respectively. 
The  tumors  can  be  seen  in  both  of  the  SART  images,  but  clutter 
shows  up  between  the  two  tumors  in  Fig.  10(e)  and  (f),  and  the 
sizes  of  the  two  tumors  in  Fig.  10(f)  are  larger  than  their  true 
sizes. 

VI.  Conclusion 

An  investigation  of  using  a  multifrequency  microwave-in¬ 
duced  TAI  system  for  early  breast  cancer  detection  has  been 
reported  in  this  paper.  The  frequency  band  for  this  system  has 
been  given  based  on  the  cutoff  frequency  of  the  human  breast. 
A  simplified  semicircular  dielectric  waveguide  mode  was  used 
to  calculate  the  cutoff  frequency  in  this  paper.  By  studying  the 
microwave  energy  absorption  properties  of  breast  tissue  and 
tumor,  we  have  shown  that  the  multifrequency  microwave-in¬ 
duced  TAI  can  offer  higher  SNR,  higher  imaging  contrast,  and 
more  effective  clutter  suppression  capability  than  the  traditional 
single-frequency  microwave-induced  TAI.  A  MART  has  been 
presented  for  image  formation.  This  data-adaptive  algorithm 
can  achieve  better  resolution  and  better  interference  rejection 
capability  than  its  data-independent  counterparts,  such  as  DAS. 
The  feasibility  of  this  multifrequency  microwave-induced  TAI 
system  as  well  as  the  performance  of  the  proposed  image  re¬ 
construction  algorithm  for  early  breast  cancer  detection  have 
been  demonstrated  by  using  2-D  numerical  electromagnetic  and 
acoustic  breast  models.  The  absorbed  microwave  energy  and  the 
thermal  acoustic  field  in  the  breast  models  have  been  simulated 
using  the  FDTD  method.  Numerical  examples  have  been  used 
to  demonstrate  the  excellent  performance  of  this  multifrequency 
technique.  More  advanced  models  are  being  developed  to  fur¬ 
ther  investigate  and  validate  the  preferential  imaging  capability 
of  the  technique  for  early  breast  cancer  detection. 
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ABSTRACT 

Delay-and-sum  (DAS)  beamforming  is  the  standard  tech¬ 
nique  for  ultrasound  imaging  applications.  Due  to  its  data 
independent  property,  DAS  may  suffer  from  poorer  resolu¬ 
tion  and  worse  interference  suppression  capability  than  the 
adaptive  standard  Capon  beamformer  (SCB).  However,  the 
performance  of  SCB  is  sensitive  to  the  errors  in  the  sample 
covariance  matrix  and  the  signal  steering  vector.  Therefore, 
robust  adaptive  beamforming  techniques  are  desirable.  In  this 
paper,  we  consider  ultrasound  imaging  via  applying  a  user 
parameter  free  robust  adaptive  beamformer,  which  uses  a 
shrinkage-based  general  linear  combination  (GLC)  algorithm 
to  obtain  an  enhanced  estimate  of  the  array  covariance  matrix. 
We  present  several  multistatic  adaptive  ultrasound  imaging 
(MAUI)  approaches  based  on  GLC  to  achieve  high  resolution 
and  good  interference  suppression  capability.  The  performance 
of  the  proposed  MAUI  approaches  is  demonstrated  via  an 
experimental  example. 

Index  Terms — Adaptive  beamforming,  Ultrasound  imaging 


I.  INTRODUCTION 

Delay-and-sum  (DAS)  beamforming  is  the  standard  tech¬ 
nique  for  ultrasound  imaging  applications.  Theoretically  this 
data  independent  approach  has  lower  resolution  and  worse  in¬ 
terference  suppression  capability  than  an  adaptive  beamformer, 
e.g.,  the  standard  Capon  beamformer  (SCB)  [1].  However,  in 
practice,  there  is  a  clear  performance  degradation  for  SCB 
when  the  covariance  matrix  is  inaccurately  estimated  due  to 
limited  data  samples  and  when  the  knowledge  of  the  steering 
vector  is  imprecise  due  to  look  direction  errors  or  imperfect 
array  calibration.  Therefore,  adaptive  beamforming  approaches 
that  are  robust  to  the  aforementioned  problems  are  desired. 

Most  of  the  early  approaches  to  robust  adaptive  beamform¬ 
ing  are  ad-hoc  techniques,  e.g.,  the  traditional  diagonal  loading 
algorithm  [2],  for  which  there  is  no  clear  way  to  choose 
the  diagonal  loading  level.  The  diagonal  loading  algorithm 
has  been  previously  applied  to  ultrasound  imaging  [3],  where 
the  diagonal  loading  level  was  set  to  be  proportional  to 
the  received  power.  The  robust  Capon  beamformer  (RCB) 
presented  in  [4],  on  the  other  hand,  can  precisely  calculate 
the  diagonal  loading  level  based  on  the  uncertainty  set  of  the 
steering  vector.  RCB  was  applied  to  ultrasound  imaging  in 
[5]  and  the  results  showed  that  RCB  can  provide  much  better 
imaging  quality  than  DAS.  However,  we  still  need  to  specify 
the  uncertainty  set  parameter  for  RCB,  which  may  be  hard  to 
do  in  practice.  To  achieve  user  parameter  free  robust  adaptive 


beamforming,  we  have  recently  devised  several  beamformers 
in  [6]  based  on  the  shrinkage  method,  which  can  compute  the 
diagonal  loading  level  automatically  without  specifying  any 
user  parameters.  Among  these  beamformers,  the  general  linear 
combination  (GLC)  algorithm  performs  well,  especially  when 
the  number  of  snapshots  is  small. 

In  this  paper,  we  present  several  user  parameter  free  ap¬ 
proaches  based  on  GLC  for  multistatic  adaptive  ultrasound 
imaging  (MAUI),  which  form  images  of  the  backscattered 
energy  for  each  focal  point  within  the  region  of  interest.  All 
the  MAUI  approaches  are  two-stage  imaging  algorithms  and 
GLC  is  employed  in  each  stage.  A  similar  idea  can  be  applied 
to  microwave  imaging  to  replace  the  user  parameter  dependent 
RCB  in  each  stage  [7].  The  complete  multistatic  data  set 
for  a  given  focal  point  can  be  represented  by  the  data  cube 
shown  in  Fig.  1.  In  one  of  the  MAUI  methods,  which  we 
refer  to  as  MAUI-1,  GLC  is  used  in  Stage  I  to  obtain  a  set 
of  backscattered  signal  estimates  at  each  time  instant.  Based 
on  these  estimates,  a  scalar  waveform  is  recovered  via  GLC 
in  Stage  II,  which  is  then  used  to  compute  the  backscattered 
energy.  An  alternative  way  of  signal  processing  in  Stage  I  is  to 
compute  a  set  of  backscattered  waveforms  for  each  transmitter, 
which  is  referred  to  as  MAUI-2.  In  addition,  we  also  consider 
a  combined  method  MAUI-C,  which  uses  the  signal  estimates 
from  both  MAUI-1  and  MAUI-2  in  Stage  I  for  the  computation 
of  backscattered  energy.  An  experimental  example  will  be 
presented  to  illustrate  the  performance  of  the  MAUI  methods. 

Notation:  The  superscript  (•)*  denotes  the  conjugate  trans¬ 
pose,  (-)T  denotes  the  transpose,  denotes  rounding  to  the 
greatest  integer  less  than  x,  E(-)  is  the  expectation  operator, 
tr(-)  is  the  trace  operator,  and  ||  •  ||  denotes  the  Frobenius  norm 
for  a  matrix  or  the  Euclidean  norm  for  a  vector.  Finally  R  >  0 
means  that  R  is  positive  semi-definite. 


MAUI-l 


MAUI-2 


Receiver  Index 


Fig.  1.  The  multistatic  data  cube.  MAUI-1  processes  the  data  set  for  a  given 
time  instant  to,  while  MAUI-2  processes  the  data  set  for  a  given  transmitter 
index  i. 
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Council  (VR). 


II.  PROBLEM  FORMULATION 

Consider  an  active  array  of  M  transducers  using  the  multi¬ 
static  (also  called  MIMO  (multi-input  multi-output)  [8])  data 
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acquisition  scheme.  Each  transducer  takes  turns  to  trans¬ 
mit  the  same  pulse  while  all  of  the  transducers  record  the 
backscattered  signals.  As  a  result,  the  received  data  set 
{P  =  t  =  0,  •••,  T— 1}  comprises  the 

A-scan  data  for  all  possible  transmitter  and  receiver  pairs  of  the 
array,  where  P ij(t)  is  the  data  sequence  of  the  backscattered 
echo  at  the  jth  transducer  due  to  transmitting  a  pulse  from  the 
ith  transducer,  and  T  is  the  number  of  data  samples  for  the 
A-scan  sequence. 

To  extend  GLC  to  the  wide -band  ultrasound  imaging  ap¬ 
plication,  we  align  the  received  signals  from  the  data  set 
{P to  each  focal  point  by  inserting  appropriate  time 
delays.  Let  r*  and  r  j  denote  the  locations  of  the  ith  transmitter 
and  jth  receiver,  respectively,  and  let  17  denote  the  location 
of  a  focal  point  in  the  imaging  region  of  interest.  The  time 
delay  due  to  the  ultrasonic  wave  propagation  from  the  ith 
transmitter  to  the  focal  point  17  and  then  back  to  the  jth 
receiver  is  approximated  as 


(3,  where 

MSE(R)  =  £{||R-R||2} 

=  \\al  -  (1  -  f3)R\\2  +  f32E{ ||R  -  R|| 2 } 

=  a2L  —  2a(l  —  (3)  tr(R) 

+(1  -  /3)2||R||2  +  (32E{ ||R  -  R||2}, 

K  =  E\y(k)y*(k)].  (4) 

The  optimal  values  for  (3  and  a  can  be  readily  obtained: 


«o  =  v(\  -  f30)  =  v  P  ,  (6) 

7  +P 

where  p  =  £{||R  -  R||2},  v  =  and  7  =  jji/I  -  R||2. 

Note  that  /3o  €  [0, 1]  and  a o  >  0. 

To  estimate  a0  and  /30  from  the  given  data,  we  need  an 
estimate  of  p,  which  can  be  calculated  as  (see  [9]  for  details): 


1 

At 


~r/ii 


(i) 


where  c  is  the  sound  propagation  speed  in  the  medium  under 
interrogation,  and  At  denotes  the  sampling  interval.  Then,  the 
time  shifted  signal  for  a  given  focal  point  of  interest  17  can 
be  represented  as 

2 Up(rf,t)  =  Pid(t  +  Tij(rf)), 

i, j  =  1,  •  •  •  M;  t  =  (),•••  ,N-  1,  (2) 

where  N  is  determined  by  the  duration  of  the  transmitted  pulse 
and  the  sampling  interval  At. 

The  problem  of  interest  here  is  to  form  an  ultrasound  image 
on  a  grid  of  points  in  the  imaging  region.  The  image  is 
formed  from  the  received  data  set  {P or  more  precisely, 
{Vi,j(rf,t)},  for  each  focal  point  of  interest. 


III.  MAUI 

The  two-stage  MAUI  algorithms  use  a  GLC-based  robust 
adaptive  beamforming  algorithm  in  each  stage.  We  first  review 
the  GLC  approach  and  then  we  show  how  to  apply  GLC  to 
the  data  set  {^(r/,  t)}  in  Stages  I  and  II  of  the  proposed 
MAUI  approaches. 


A.  GLC 

In  the  GLC  approach,  we  replace  the  sample  covariance  ma¬ 
trix  in  SCB  by  an  enhanced  estimate  obtained  via  a  shrinkage- 
based  method.  The  enhanced  covariance  matrix  estimate  R  is 
obtained  by  linearly  combining  the  sample  covariance  matrix 
R  and  a  shrinkage  target  (we  use  the  identity  matrix  I  here 
for  lack  of  a  better  choice)  in  an  optimal  mean-squared  error 
(MSE)  sense: 

R  =  al  +  (3  R,  (3) 

where  R  =  J2k=i  y(^)y*(&)>  with  the  L  x  1  vector  y(k) 
denoting  the  kth  snapshot  and  K  representing  the  total  number 
of  snapshots.  The  shrinkage  parameters  a  and  /3  in  (3)  are 
estimated  by  minimizing  the  MSE  of  R  with  respect  to  a  and 


P 


1 


f>(fc)||4-E||R 


k= 1 


(7) 


Using  (7)  we  can  get  estimates  for  /30  and  a0  as 


A)  — 


7 


and 


7  +  P 
q/q  =  0(1  -  (30), 


(8) 

(9) 


where  0  =  tr^ ,  and  7  =  || z>I  —  R||2.  Note  that  do  >  0  and 
(3q  >  0,  which  guarantees  that  the  enhanced  covariance  matrix 
estimate  R  >  0. 

Substituting  (8)-(9)  in  (3)  yields  the  enhanced  covariance 
matrix  estimate  R.  Using  R  instead  of  R  in  the  SCB  formu¬ 
lation,  we  obtain  the  beamformer  weight  vector  for  GLC  as 
follows: 


w 


R  a 

a*R_1a 


(10) 


where  a  denotes  the  assumed  steering  vector  [10].  Note  that 
GLC  is  a  diagonal  loading  approach  with  the  diagonal  loading 
level  a0//3o  determined  automatically  from  the  observed  data 
snapshots  {y (k)}£=1. 


B.  Stage  I 

To  apply  the  GLC-based  robust  adaptive  beamformer  to  the 
data  set  {^y(r/,t)}  in  (2),  we  use  two  approximate  signal 
models  for  yi,j(rf,t)  by  making  different  assumptions.  Since 
we  will  concentrate  on  the  focal  point  17  in  what  follows,  the 
dependence  on  7  will  be  dropped  for  notational  simplicity. 

The  MAUI-1  algorithm  uses  the  following  signal  model: 

y  i(t)  =  a  (t)si(t)  +  ei(t ),  (11) 

where  y i(t)  =  {2/^,1  (t) ,  •  •  •  ,^,m(7)]T  represents  the  aligned 
array  data  vector  of  the  ith  transmitter,  Si(t)  denotes  the 
signal  of  interest  (SOI)  that  is  proportional  to  the  ultrasound 
reflectivity  or  scattering  strength,  which  is  assumed  to  depend 
on  the  transmitter  i  but  not  on  the  receiver  j,  ei(t)  denotes  the 
residual  term  due  to  noise  and  interferences,  and  a (t)  denotes 
the  array  steering  vector  that  is  assumed  to  be  approximately 
equal  to  Imxi-  Here,  we  assume  that  a (t)  may  vary  with  t, 
but  is  constant  with  respect  to  the  transmitter  index  i. 
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In  Stage  I,  for  a  given  time  to,  we  form  a  pseudo-covariance 
matrix  by  considering  the  number  of  transmitters  as  the 
number  of  snapshots: 

RVo)  =  LY(f0)Y*(f0), 

Y(t0)  =  [yi(t0)---yM(t0)].  (12) 

By  using  R(t0)  as  the  sample  covariance  matrix  we  obtain 
an  enhanced  covariance  matrix  estimate  R (t0)  as  described  in 
Section  III. A,  and  then  calculate  the  weight  vector  w (to)  for 
Stage  I  of  MAUI- 1  using  (10)  with  a  =  lMxi  as: 


w (t0)  = 


R(tp)  *a 
a*R(t0)_1a 


(13) 


Once  we  got  the  weight  vector,  we  can  estimate  Si(t0)  in  (11) 
as: 


Si(t0)  =  w*(t0)y;(to).  (14) 


Define  a  vector  s(t0)  =  [si(^o)5'"  ,%(to)]T  of  the  esti¬ 
mated  signals  for  all  transmitters.  Repeating  the  above  process 
from  t0  =  0  to  t0  =  tV  —  1,  we  build  the  matrix  Si  = 
[s(0),  •  •  •  ,s(N  —  1)]. 

The  MAUI-2  algorithm  considers  another  signal  model: 


y  i(t)  =  a  +  ei(t),  (15) 


where  a^  denotes  the  array  steering  vector,  which  is  also  as¬ 
sumed  to  be  approximately  equal  to  Imxi-  However,  different 
from  MAUI-1,  here  a i  is  assumed  to  change  with  i,  but  be 
constant  with  respect  to  t. 

For  a  given  transmitter  i,  the  covariance  matrix  in  Stage  I 
of  MAUI-2  is  formulated  as: 


_ YY* 

N  %  ** 

[y-i(o)  -  -  -  yi(JV  - 1)] . 


(16) 


Using  Rj  as  the  sample  covariance  matrix  we  get  an 
enhanced  estimate  R$,  and  then  compute  a  weight  vector 
w i  using  (10).  The  time  sample  vector  of  the  corresponding 
beamformer  output  can  be  written  as 

Si  =  [w*Y  if.  (17) 


C.  Stage  II 

In  Stage  II,  the  signal  model  for  both  MAUI-1  and  MAUI-2 
can  be  represented  as: 

^m(t)  &mS(j2)  T  ®ra(^)j  t  —  0,  •  •  •  ,  1,  171  =  1,  2, 

(18) 

where  the  steering  vector  a m  is  assumed  to  be  1mxi>  and 
em(t)  represents  the  residual  term.  Similar  to  Stage  I,  the 
knowledge  of  a m  may  be  imprecise  and  the  sample  size 
N  may  be  small.  Hence  the  GLC-based  robust  adaptive 
beamformer  is  used  again  to  estimate  s(t).  Taking  Rm  as  the 
sample  covariance  matrix: 

1  7V_1 

*lm=  N  12  m  =  1,2,  (19) 

t= 0 

and  paralleling  the  development  in  Stage  I,  we  obtain  the 
weight  vector  wm  using  (10).  Then,  the  output  signal  estimate 
is  computed  as: 

s(t)  =  W*msm(t),  m=  1,2.  (20) 

Finally,  the  signal  energy  for  a  particular  focal  point  17  is 
computed  as: 

N-l 

£(rf)  =  ^2  s2(t).  (21) 

£=0 

For  Stage  II  of  MAUI-C,  the  signal  model  can  be  written 
as: 


s c(t)  =  a cs(t)  +  ec(t),  t  =  0,  •  •  •  ,  N  -  1,  (22) 

where  the  vector  s c(t)  is  considered  now  to  be  a  snapshot 
from  a  2M-element  “array”,  and  the  steering  vector  a c  is 
assumed  to  be  12mxi-  &c{t)  denotes  the  residual  term.  We 
obtain  the  weight  vector  w c  for  MAUI-C  via  (10)  by  using 
the  following  sample  covariance  matrix: 

1  Ar_1 

Rc  =  jz  22  sc(t)s*c(t).  (23) 

t= 0 

The  beamformer  w c  yields  an  estimate  of  the  signal: 


Repeating  the  above  process  for  ,  M  yields  a  set  of 

waveforms  S2  P  [si,  •  •  •  ,  s m]T- 

As  we  mentioned  before,  the  errors  in  the  sample  covariance 
matrix  and  the  steering  vector  cause  performance  degradations 
for  any  adaptive  beamforming  algorithms.  GLC  is  designed  to 
improve  the  covariance  matrix  estimate.  MAUI-1  and  MAUI-2 
use  different  sample  covariance  matrices.  Hence  the  improve¬ 
ments  obtained  by  using  GLC  may  be  different.  This  fact 
motivates  us  to  combine  MAUI-1  and  MAUI-2  to  achieve  a 
better  performance.  We  refer  to  this  combined  method,  where 
Si  of  MAUI-1  and  S2  of  MAUI-2  are  used  simultaneously, 
as  MAUI-C.  We  denote  the  combined  signal  matrix  as  Sc  = 

[s T  *f. 

Let  the  M  x  1  vectors  {sm(t)}t=o,—  ,n-i  denote  the 
columns  of  Sm  for  m  =  1,2,  and  let  the  2 M  x  1  vectors 
{sc(t)}t=o,-  ,n-  1  denote  the  columns  of  S c-  Not e  that  both 
MAUI-1  and  MAUI-2  obtain  M  signal  waveform  estimates 
at  the  end  of  Stage  I,  while  MAUI-C  obtains  2 M  signal 
waveform  estimates.  We  will  apply  GLC  to  these  estimates 
in  Stage  II  to  recover  a  scalar  waveform  and  compute  the 
signal  energy  at  the  focal  point. 


s(t)  =  wj.s  c(t).  (24) 

Then,  the  backscattered  energy  at  the  focal  point  17  is  com¬ 
puted  via  (21). 

IV.  EXPERIMENTAL  EXAMPLE 

In  this  section,  we  present  some  experimental  results  to 
demonstrate  the  performance  of  the  three  MAUI  algorithms. 
The  complete  multistatic  data  set  was  obtained  by  Bioacoustics 
Research  Laboratory  of  the  University  of  Illinois  at  Urbana- 
Champaign.  The  scene  of  interest  contains  several  wire  targets 
arranged  in  a  complicated  pattern.  The  data  was  collected  us¬ 
ing  a  64-element  linear  array.  The  transducer  center  frequency 
was  2.6  MHz,  the  sampling  rate  was  25  MHz,  and  the  sound 
velocity  was  assumed  to  be  1450  m/s.  For  comparison,  the 
multistatic  DAS  scheme  is  also  applied  to  the  same  data  set. 
The  DAS  scheme  estimates  the  signal  waveform  s(t)  as 

s(t)  =  w*asY (t) wDAS,  t  =  0,  •  •  •  ,  N  -  1,  (25) 

where  wDAS  =  a/M  is  the  weight  vector  for  DAS.  The 
backscattered  energy  at  17  is  then  estimated  via  (21). 
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Fig.  2  shows  the  ultrasound  images  for  the  wire  data  set 
under  consideration.  The  images  are  displayed  on  a  logarith¬ 
mic  scale  with  a  30  dB  dynamic  range.  In  Figs.  2  (a)-(d),  we 
compare  the  images  obtained  via  DAS  and  MAUI  algorithms 
using  only  the  central  32  elements  of  the  array  (M  =  32). 
Since  DAS  simply  sums  all  signals,  the  DAS  image  shown 
in  Fig.  2  (a)  has  higher  sidelobe  level  and  poorer  resolution 
than  the  MAUI  images  shown  in  Figs.  2  (b)-(d).  Comparing 
Fig.  2  (b)  and  Fig.  2  (c),  which  correspond  to  MAUI-1  and 
MAUI-2  respectively,  we  note  that  MAUI-2  image  has  a 
lower  background  clutter  level.  However,  MAUI-2  has  poorer 
resolution:  some  wire  targets  are  not  discernable  in  the  MAUI- 
2  image.  On  the  other  hand,  the  image  obtained  via  MAUI-C 
has  low  sidelobe  level  similarly  to  MAUI-2  and  high  resolution 
similarly  to  MAUI-1.  Moreover,  all  targets  are  clearly  shown  in 
the  MAUI-C  image.  For  comparison,  we  also  include  the  DAS 
image  obtained  using  the  entire  array  (M  =  64).  Note  that 
MAUI  algorithms,  especially  MAUI-C,  with  32  transducers 
can  achieve  similar  imaging  quality  to  DAS  with  a  double 
sized  array. 


V.  CONCLUSIONS 

We  have  presented  three  user  parameter  free  approaches  to 
multistatic  adaptive  ultrasound  imaging  (MAUI).  These  two- 
stage  MAUI  approaches  employ  a  GLC-based  robust  adaptive 
beamformer  in  each  stage  to  achieve  high  resolution  and  good 
interference  suppression  capability,  and  also  they  are  robust  to 
small  sample  size  problems  and  array  steering  vector  errors. 
More  importantly,  GLC  is  a  user  parameter  free  approach 
as  opposed  to  other  existing  robust  adaptive  beamforming 
algorithms,  which  makes  it  easy  to  use  it  in  practice.  The 
experimental  results  have  demonstrated  the  effectiveness  of  the 
MAUI  algorithms  for  ultrasound  imaging.  We  have  shown  that 
the  MAUI-C  method,  which  combines  MAUI-1  and  MAUI-2, 
provides  the  best  imaging  quality. 
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Fig.  2.  Ultrasound  images  obtained  via  (a)  DAS  with  M  =  32,  (b)  MAUI-1 
with  M  =  32,  (c)  MAUI-2  with  M  =  32,  (d)  MAUI-C  with  M  =  32,  and 
(e)  DAS  with  M  =  64. 
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